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Abstract 

The modifiable areal unit problem and the ecological fallacy are known problems that occur when 
modeling multiscale spatial processes. We investigate how these forms of spatial aggregation er¬ 
ror can guide a regionalization over a spatial domain of interest. By “regionalization” we mean a 
specification of geographies that define the spatial support for areal data. This topic has been stud¬ 
ied vigorously by geographers, but has been given less attention by spatial statisticians. Thus, we 
propose a criterion for spatial aggregation error (CAGE), which we minimize to obtain an optimal 
regionalization. To define CAGE we draw a connection between spatial aggregation error and a 
new multiscale representation of the Karhunen-Eoeve (K-E) expansion. This relationship between 
CAGE and the multiscale K-E expansion leads to illuminating theoretical developments includ¬ 
ing: connections between spatial aggregation error, squared prediction error, spatial variance, and 
a novel extension of Obled-Creutin eigenfunctions. The effectiveness of our approach is demon¬ 
strated through an analysis of two datasets, one using the American Community Survey and one 
related to environmental ocean winds. 
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1 Introduction 


There has long be en interest i n non- st atistieal method s 


spatial data (e.g., 


Openshawl (Il977h . 


MurtaghI (Il992h . 


or spee i 


Martini (l2002h . 


ymg 


g eographies t o sum marize 


Guol (l2008h . and 


Logan 


(I 2 OI ll) l. In general, this is known as “regionalization,” and it is an important (and sometimes 


required) task for many applieations. For example, the Ameriean Community Survey (ACS) is 
an ongoing survey administered by the US Census Bureau that produces estimates of important 
US demographic variables. The ACS provides public-use data referenced over areal units (e.g., 
median household income over US counties). Similar to the decennial census, many of these ge¬ 
ographic regions are required (e.g., states, counties, etc.), however, other regions are consistently 
being evaluated and changed (e.g., combined statistical areas, metropolitan divisions , metropoli 


tan sta tistical areas, etc.) in a sub-optimal manner based on population controls (e.g.. 


Blank et al. 


(I 2 OI lll l. This suggests that there is a clear need for regionalization methodology. Thus, we develop 


a principled statistical methodology for evaluating spatial aggregation error and optimal statistical 
regionalization. 

Regionalization is a topic that has been considere d primarily by geographers. The current 
state-of-the-art is the deterministic “max-p algorithm” (Duoue et ah. 


2013 


Folch and Spielman. 


2014; 


Spielman and Logan. 


2012 


2015b . In general, the max-p algorithm is a 


greedy search algorithm (using any desired criterion) that groups data defined on areal units into 
P (< ^a) contiguous regions. The max-p algorithm offers a solution, but there are many known 
pitfalls to this approach. The most significant issue from the perspective of multiscale spatial in- 
ference is that th e regions obtained by this approach are not protected from the ecological fallacy 


Robinson. 


1950b . Hence, proper inferential conclusions must be limited to a single (often difficult 


to interpret) spatial support. 

We interpret the ecological fallacy as a type of spatial aggregation error, which will be criti¬ 
cal to our approach for regionalization. In particular, the ecological fallacy refers to the situation 
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where eonelusions at the point-level spatial support differ from eonclusions at an aggregate-level 
spatial support. Similarly, ecological inference is explicitly defined as inference on individual be¬ 
havior drawn from aggregate data (also sometimes referred to as downscaling). This topi c has 


experi enced growing interest within a v ariety of subiect ma tter disciplines. For example, see 


(1 19971) for the sociological d ata setting; 


cations in epidemiology; and 


D^bv et al. 


Mearnsetal 


King 


(l200ll) . and the references therein, for appli- 


(12014). and the referenc es therein, for the climatology 


setting. Following the terminology of 
as image segmentation , whic h involves optima l 


Kolaczvk and Huangl (1200 lb . a similar problem is known 


Kolaczvk and NowakI (12004) . 


Kolaczvk et al 


V div iding an image into smaller regions (e.g., see 


eco 

and 


ogical inference and im age segmentation see 


(l2005b. and 


Ferreira et al 


Ferreira and Led ( 2007 ). 


Wakefieldl (12004) . 


(1201 lb). For rev i ews o f 


Waller and Gotwavl (12004) . 


The mo diliahle areal unit nrohlem tMAUP) is another tvne of spatial aggregation error. Waller 


and Gotway (12004) consider the MAUP to be the geographic manifestation of the ecological fal¬ 


lacy. That is, the MAUP refers to situations where conclusions on one aggregate spatial support 
differ from conclusions on another distinct aggregate spatial support, whereas, the ecological fal¬ 
lacy concerns conflicting conclusions at point-level and aggreg ate-level supp orts. The MAUP has a 


rich hi story, originally considered by 


1979 ). 


199.3b . 


Gehike and Biehll (119.34) . and later by 


Openshaw and Tavloi 


Recently, the MAU P has become a topic covered in s tanda rd textbooks inc 


Waller and Gotwayl (12004) . 


Cressie and Wikld (1201 lb . and 


Banerjee et al. 


uding 


Cressie 


(1201.5b . among 


others. 

The aforementioned forms of spatial aggregation error are closely related to the spatial change 


of support (COS) problem, which refers to conducting statistica 


fro m the spatial support of the data (e.g.. 


and 


Banerjee et al. 


Waller and Gotwavl (12004) . 


infer e nce on a support tha t differ s 


Cressie and Wikld (1201 lb . 


(12015b ). Methods for spatial COS allow one to choose any support on which to 


perform statistical inference. However, different choices for the spatial support result in different 
magnitudes of spatial aggregation error. Nevertheless, the inherent flexibility to use any desired 
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spatial support for inference has made spatial COS a popular area of research in both mu l tiscale 


spatial analysis and other subject m atter disciplines . For e xample, see 


for t he environmental data setting: 


Mugglin et al. 


Wikle and Berlinen (120051) 


(1 1998b for the public health setting: Bradley 


et al. (2015bb for the survey data setting; and 


Waller and Gotwavl (12004) and 


Trevisani and Gelfand 


(120 13b for a review. To capitalize on the flexibility of spatial COS methods, we adopt a multiscale 


spatial perspective to quantify spatial aggregation error and to develop a method for regionaliza¬ 
tion. 

The known presence of spatial aggregation error suggests an approach for an optimal region¬ 
alization. Specifically, our primary inferential question is the following: can we choose a spatial 
support that minimizes spatial aggregation error? To motivate this perspective, consider an exam¬ 
ple dataset obtained from the ACS. In Figures 1(a) and 1(b), we plot 5-year period estimates of 
median household income by county and state, respectively, for 2013. Upon comparison. Figures 
1(a) and 1(b) show that the state-level ACS estimates suffer from noticeable spatial aggregation 
error. For example. Figure 1(b) suggests that households in Virginia have moderately high in¬ 
come, yet Figure 1(a) shows that only households in counties near Richmond have high income. 
Similarly, Figure 1(b) suggests that households in New York state have a moderately high income 
while Figure 1(a) shows that only households in counties near Manhattan have high income. These 
examples, and many others that are quite obvious upon study of these figures, provide evidence 
that states are not an appropriate (i.e., optimal) spatial support to summarize median household 
income, political reasons notwithstanding. 

In what follows, we formalize this intuition and develop a criterion to quantify spatial aggre¬ 
gation error and an associated method for regionalization. Our approach is to quantify spatial 
aggregation error using what we call the criterion for spatial aggregation error (CAGE). Hence, 
an optimal spatial support is obtained by minimizing CAGE. T he primary theoretical tool use d 


to develop this criterion is the Karhunen-Eoeve (K-E) expansion (IKarhunen . 


1947 


Eoeve. 


1978b . 


which is a well-known representation of a point-referenced process as the weighted sum of spa- 


3 



































Figure 1: ACS 5-year period estimates of median household ineome for 2013. In (a), we plot the 
ACS estimates by counties, and in (b) we plot the ACS estimates by state. We superimpose the 
state boundaries as a reference in both panels. Notice that the color-scales are different for each 
panel. In (b), the borders of the states are highlighted in white except for New York and Virginia, 
whose borders are highlighted in black. Also, Richmond Virginia and Manhattan are indicated 
with arrows in (b). 


tially varying eigenfunctions, where the weights are random. In more precise terms, we develop 
CAGE through a powerful technical result, which dictates that spatial aggregation error does not 
occur when the eigenfunctions of a spatial random process are constant between spatial scales. 
Thus, CAGE is a measure of between spatial scale homogeneity of eigenfunctions within a novel 
multiscale representation of the K-E expansion. 

To date, there has been no such criterion that quantifies spatial aggregation e rror in this ma n- 
ner. The spatial statistics literature places an emphasis on prediction error (e.g., ICressid (119931) ). 


and thus, such an aggregation-based approach for uncertainty quantification offers an exciting new 
perspective for spatial statistics. Therefore, to develop this perspective we provide technical results 
relating CAGE to prediction error and spatial variance. 
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After having defined CAGE, we can choose a regionalization in a manner that mitigates spa¬ 
tial aggregation error. In particular, we propose an efficient search algorithm (with CAGE as the 
selection criterion) to specify a regionalization over the spatial domain of interest. This search 
algorithm involves two stages. In the first stage, a naive algorithm, sav ^-means (e.g.. Hartigan 


and Wong ( 


1979t) ) is used to determine a collection of spatial supports from which to select. Then, 


in the second stage CAGE is used to select a single spatial support from among the collection of 
spatial supports determined in the first stage of the search algorithm. This two-stage approach is 
extremely efficient because it uses an easy-to-compute deterministic algorithm to direct the path of 
spatial supports from which to choose. As such, it can be incorporated efficiently within a Bayesian 
framework using a Markov chain Monte Carlo (MCMC) implementation of a latent spatial model, 
which facilitates uncertainty quantification. 

Einally, to apply our search algorithm in practice, we provide a specification for the multiscale 
eigenfunctions. Thus, we introduce a general class of eigenfunctions that leads to a consistent 


class of multi scale spatial processes. T o do this, we utilize t 


framework of 


Obled and CreutinI (Il986h . 


le often overlooked, but remarkable 


Obled and CreutinI (Il986h show that any class of geosta- 


tistical basis functions can be re-weighted so that they are eigenfunctions within a (single-scaled) 
K-E expansion. This notion of what we call generating basis functions (GBEs), is central to our 
development of multiscale eigenfunctions. As interest in spatial and spatio-temporal processes has 
turned to “big data” problems with large numbers of prediction and/or data locations, the model¬ 
ing focus has shifted to th is basis function persp ective incorporating complete, over-complete, and 


reduced-rank expansions (IBradlev et ah . 


2015ah . Thus, the use of GBEs greatly increases the gen¬ 


erality and utility of our approach. Eurthermore, the use of GBEs is a necessity for our approach 
to regionalization because they allow us to perform spatial COS without assuming some form of 
between scale homogeneity. 

The remainder of this paper is organized as follows. In Section 2, we introduce the multiscale 
K-E expansion and CAGE. Next, in Section 3 we describe how to use CAGE in practice, which 
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includes details on truncating the multiscale K-L expansion and the introduetion of the two stage 
regionalization algorithm. Seetion 4 provides derivations of a eonsistent elass of multiscale eigen- 
funetions to use within the CAGE framework. Then, in Seetion 5 a demonstration is given using the 
motivating dataset of ACS 5-year period estimates of median household ineome from Figure 1. In 
addition to demonstrating the regionalization algorithm for ACS period estimates, this applieation 
also highlights an important use of optimal regionalization, namely, aggregation for the purpose of 
dimension reduction. Finally, Section 6 contains a concluding discussion. We provide additional 
Supplemental Materials including: the proofs of teehnieal results, simulation studies, and an addi¬ 
tional applie ation using a datase t eonsisting of Mediterranean wind measurements (a subset of the 


data used in 


Milliff et al. 


(1201 ih l. The Mediterranean wind example is used to illustrate that the 


two-stage regionalization algorithm is flexible enough to handle multiseale spatial data. 


2 Quantifying Aggregation Error 

Here, we provide requisite extensions of the K-F expansion to the multiscale setting (Section 2.1). 
These results are then used to formally define CAGE (Seetion 2.2). 

2.1 The Multiscale Karhunen-Loeve Expansion 

Consider a real-valued spatial proeess that is realized at (possibly) both point-level and aggregate- 
level spatial supports. That is, the values in the sets {Tj(s) : s G and {Ta(A) : A G Da} can 
be realized, where is a eontinuous spatial random process defined on C and Ya is 

a spatial random proeess defined on areal support Da with Da = {A, : i = 1,..., and A, C 
The set A, is an areal unit (e.g., a eounty, state, or census tract) and may be overlapping, contained 
in, or superimposed over another distinet areal unit A/ G Da for j ^ i. 
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The corresponding multiscale spatial process can be written as 


y(u) = 


K(u) ifueD, 


}a(u) ifue Da', ueDsUDa. 


( 1 ) 


We interpret iA(') as being com puted from the point-level process particular, as is 

standard in spatial statistics (e.g., 


Cressid (119931) . p. 284), we assume 


YAiA) = -^ [ y,(s)Js; AeDa, 

\A\ Ja 


( 2 ) 


where |A| represents the cardinality of the set A. Consequently, placing a statistical model on y^ 
implicitly places a statistical model on Ya and Y through ([U) a nd (l2l). We explore t his de pendency 


between ([T]) and dH) using the well-known K-L expansion (e.g.. 


Cressie and Wikld (1201 ll) . p. 156), 


= L0,/(s)a;; SGD„ 

f=l 


(3) 


where, without loss of generality, {y^(-)} is assumed to be mean-zero, the random variables in the 
set {«/: 7 = 1,2,...} are uncorrelated with associated variances {A/: y = 1,2,...} (called eigen¬ 
values), the orthonormal real-valued functions {^j{s) : j = 1,2,...} (called eigenfunctions) have 
domain D^, and satisfy a Fredholm integral equation for a given valid covariance function. (Note 
that the conditions needed for the K-L expansion are given in the statement of Proposition 1.) 

The use of the K-L expansion greatly increases the generality of our approach, since Mercer’s 


theorem dictates t 


expansion (|Mercei , 


l at poi nt-level covariance functions can be dec omposed according to the K- 


1909h under a very general set of assumptions (|Ferreira and Menegatto . 


2009 


This leads us to define a multiscale K-L expansion, which we formalize through Proposition 1 
below. 
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Proposition 1: Let S^') be a probability space, where Q is a sample space, ^ is a sigma- 

algebra on Q,, and LP is a finite Borel measure. Let Ys{s) be defined by the mapping : Ds x Q. 

M, such that 7^(5) is measurable for every s G D^, and Dg C is a topological Hausdorff space. 
Assume thatC{s,u) = cov { 7 ^( 5 ), y^(H)} is a valid covariance function that exists for each s,u G 
Let L^(Q) denote the Hilbert space of real-valued square integrable random variables. 

i. Then, for each A G Dgwe have that 

00 

YA{A) = Y,^A,M)ap (4) 

i=i 

in L?{Q.), where for each positive integer j, ^A,j{A) = JA^j{s)ds/\A\, the random variables 
in the set {aj : 7 = 1 , 2 ,...} are uncorrelated with associated variances {Xj : j = 1 , 2 ,...} 
(called eigenvalues), the orthonormal real-valued functions : j = 1 , 2 ,...} (called 

eigenfunctions) have domain D^, and satisfy the Fredholm integral equation for C{s,u). 

a. Then for any A G Ds and B G Ds we have that 

cov{7a(A),7a( 5)}= lim (5) 

i=\ 

The proof of this proposition can be found in the Supplemental Materials. 

Remark 1: We call the expression in ® the multiscale K-L expansion since Proposition l.i ex¬ 
tends the K-L expansion in ® to a similar infinite-dimensional process that is a function of any 
A G Ds. Similarly, the expression in ® can be seen as an extension of Mercer’s theorem to the 
multiscale spatial setting. 

Remark 2: In practice, the latent multiscale spatial process of interest 7 is not observed perfectly. 
Instead, we observe the n-dimensional data vector given by Z = (Z(u) : u G D^UD^)', where the 


8 


observed locations are denoted by = {sf : i = 1, } C Ds and = {Aj : j = 1, C 

Da, and n = n^ + n^. We assume that the stochastic processes Z : x —)■ M and Y are generated 
based on the generic probability space ^such that the conditional probability density 
function of Z(u) |Z exists for each u G UD^. 

Remark 3: For purposes of implementation it is helpful to define a set Dg = {Bj : j = 1, ...,ns} 
with Bj = 0 for j ^ and Bj C Ds for each j. Here, Dg represents the finest resolution 
spatial support on which one is willing to perform inference. Then, after observing data Z(-), 
statistical inference is performed using sample draws from the distribution of Y^jZ, where the 
ng-dimensional process vector is given by Yg = {Ya{B) : B G Dg)'. 

2.2 The Criterion for Spatial Aggregation Error (CAGE) 

There is an implicit conceptual challenge involved with quantifying spatial aggregation error. As 
Gotway and Waller (2011) discuss, the consequences of spatial aggregation error extend beyond 
between-scale differences of the values of a single statistic (e.g., correlation coefficient, mean, 
etc.). Thus, we say that spatial aggregation error occurs when there are between-scale differences 
for any generic statistic. The multiscale K-L expansion in ® provides insight on a formalization 
of this concept, which we state in Proposition 2. 

Proposition 2: Assume that the conditions of Proposition 1 hold. Let f be a measurable real¬ 
valued function with domain that is discontinuous only on a set with measure zero. Let Xk 
be strictly greater than zero for each k = 1,2, .... Define a generic point-level support {xj : j = 
l,...,nA}, such thatXj G Bj C Aj G Da for j = \,...,nA, = [Ys{xj) : j = \,...,n/f)', Y^^^ = 

(YA{Bj) : j = 1,...,^^)^, andYA = (Fa(A) :A ^Da)'. Then the following statements hold for Y {■) 
in ([T]).- 
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i. ^kixj) = for j = 1, ond every positive integer k, if and only if = fiY^) 

almost surely. 

a. (j)k{Bj) = (j)k{Aj)for j = 1, ond every positive integer k, if and only if f{Y^g^) = fiY^) 
almost surely. 

Hi. If^kftj) = (j)k{Aj) for every positive integer k, and every xj E Bj and j, then f{Y^^'’) = f{YA) 
almost surely. 

(A) 

Remark 4: Proposition 2 provides a condition so that there is no ecological fallacy between ’ 
and Y^, and no MAUP between Y^^ and Y^. By “no ecological fallacy” and “no MAUP,” we 
mean that for any real-valued, measurable, (almost) continuous statistic /, = /(Y^) and 

/(Y^^) = /(Y^) almost surely. This ensures that conclusions using the summary statistic / stay 
the same regardless of the scale of Y. In general terms. Propositions 2.z and l.ii show that “no 
spatial aggregation error” is equivalent to between-scale homogeneity of eigenfunctions within a 
multiscale K-L expansion. Furthermore, Propositions l.i and l.iii provide a relationship between 
the ecological fallacy and the MAUP; namely, if there is uniformly no ecological fallacy for any of 
the sets in {Bj} (i.e., = ^{Aj) for every xj E Bj and j), then there is no MAUP. 


Proposition 2 guarantees that spatial aggregation error does not occur when the point-level 
eigenfunctions are constant over each region in Da. This leads naturally to a criterion that measures 
departures from the absence of spatial aggregation error. Specifically, we define CAGE as follows: 


CAGE(A) = E 


\A\ 


■ds\X 


( 6 ) 


where A is a generic areal unit (i.e., A C Ds), and the expectation is taken with respect to the con¬ 
ditional distribution given the data. The logic behind © is straightforward: if CAGE(A) is equal 
to zero there is no loss of information when aggregating to Da, and if CAGE(A) is close to (far 
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from) zero then we lose a small (large) amount of point-level information when aggregating to A. 
Hence, maps of {CAGE(A/) : i = 1, can be used to assess whether statistical inference on 

Ya is reasonable relative to the point level process. 


In some settings the latent process cannot realistically be defined at the point level. For ex- 


ample, the mec 


Banerjee et al. 


i an (ov er counties) household income in Figure 1 cannot be interpreted on Ds (see 


(12015b for a discussion and more examples). Hence, for these settings the mul¬ 


tiscale K-F expansion is used for spatial change of support, and the lowest spatial resolution on 
which Y is defined is Dg. We use the following discretized CAGF (abbreviated as “DCAGE”) in 
these settings: 


DCAGE (C) = E 


heH 


C 


(7) 


where C = // C {1, ...,^ 5 }, and G Dg for each h e H. Proposition 2.ii implies the 

following logic for ([7]): if DCAGE (C) is equal to zero there is no loss of information when ag¬ 
gregating Dg to higher spatial resolutions, and if DCAGE(C) is close to (far from) zero then we 
lose a small (large) amount of lower resolution information when aggregating Dg to higher spatial 
resolutions (see Remark 3). 

To date there has been no attempt to quantify the magnitude of spatial aggregation error using 
criteria like db]) and ([7]). In the geostatistic al setting, emphasis is usually placed on minimizing 


the squared prediction error dCressie . 


1993b . From this point-of-view, it is worthwhile to note that 


there are connections between the squared prediction error, spatial variance, and CAGE in d^, 
which we formally state in Proposition 3 below. 


Proposition 3: Assume that the conditions of Proposition 1 hold. Also, assume that the stochastic 
process Z : x —)■ M A generated based on a generic probability space (Q, 3^, 0^') such that the 
conditional probability density function ofY {u)\Z exists for each u G UD^, where Z is defined 
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in Remark 2. Then, CAGE in ® has the following alternative expressions: 


CAGE {A) =E 
CAGE {A) =E 

CAGE{A)=E 


Ja W 




{n(s)-?4(A))" 


( ^ '1 2 

^ ds:\7 

|A| 

-E 

|Ta(A)-Ta(A)| |Z 


( 8 ) 

(9) 

( 10 ) 


where A is a generic areal unit (i.e., A C Ds), and Ya (A) = E[Ya (A) |Z). 


Remark 5: Each expression in Proposition 3 provides interesting motivation for CAGE. Eor ex¬ 
ample, © was motivated by Proposition 2 (i.e., by measuring the departure from the absenee of 
spatial aggregation error), however, one eould argue to use ® from a praetieal perspeetive. That 
is, intuition suggests that it is reasonable to make finer seale inference using the aggregate proeess 
if Ef (s) is eonsistently “elose” to Ya (A). However, it is important to note that our use of the K-E ex¬ 
pansion is important beeause it allows us to perform spatial ehange of support to obtain Ya without 
assumptions of between-seale homogeneity. Additionally, the expression in ® is especially inter¬ 


esting from a historical perspective , since many of t 


foeused on second order statisties (IRobinson . 


le early referenees on spatial aggregation error 


195flt) . Here, we see that between-scale differences 


of varianees have a eonneetion (through Propositions 1, 2, and 3) to between-seale differenees of 
any statistic. 


Remark 6: The “ANOVA-type” decomposition in ® offers a different perspective in which to 
interpret ®. The first term on the right-hand-side of ® (from left to right) represents a within- 
areal unit predietion error. Speeifioally, the first term represents the predietion error between the 
point-level proeess Ys and the aggregate-level estimator Ya- The seeond term in ® shows that a 
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minimax-type approach is used for between areal unit error. That is, we minimize the squared 
prediction error to obtain Ya, but penalize for choosing A so that Ya is close to Ya. One could 
conceive of a version of Proposition 3 that provides similar identities for the DC AGE in (|7]). In 
Supplemental Materials, we provide the statement and proof of this technical result. 


3 Statistical Methodology for Regionalization 


In practice, higher order components, of the infinite sum in ® , correspond to a decreasing percent¬ 
age of variation. Thus, it is st andard practice to truncate the K-L expansion, and assum e that the 


residual is negligible (e.g., see 


Obled and CreutinI (Il986h . and 


Cressie and Wikld (1201 ih p. 267). 


In this section, we extend the results from Section 2 to accommodate this common assumption. 
In particular, for our applications we truncate the multiscale K-L expansion (Section 3.1), which 
leads to another version of CAGE (Section 3.2). With these details in place, we can describe how 
to use CAGE for regionalization (Section 3.3). 


3.1 The Truncated Multiscale Karhunen-Loeve Expansion 

A common simplification of the K-L expansion is to truncate the infinite sum in ® and assume 
that 


<l>s) = L 8 e D„ (11) 

y=i 

where r is a fixed and “known” integer, the r-dimensional vector of eigenfunctions is given by 
0 ^(-) = ( 0^,1 (•), and the associated r-dimensional random vector is a = (tti,..., Ct^)'. 

It is important to note that Tj(s; 0^) ^ Tj(s) in general due to the truncation in (fTTI) . 
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Now, © and (fTTl) provide an immediate expression for Ya, namely, 



aj = <l>{A;<l>^ya; AeDa, 


( 12 ) 


where ^(A; 0^) = Jj^^s,j{s)ds : j = l,...,rj . Then, ©, (fTTl) . and (IT^ imply the following 
expression for the truneated K-L expansion of the multiseale spatial proeess Y, 




^(u;0J'a if uE Da; seDsUDa, 


if u G 


(13) 


where it is important to note that the r-dimensional random veetor a is the same for both supports. 
Validity of the implied eovarianee funetion for Y follows immediately from the quadratie form (see 
Supplemental Materials for more details). 

The distributional assumptions governing Propositions 1—3 were very general (see Remark 2). 
For the truneated multiseale K-L expansion we ineorporate additional distributional assumptions. 
In partieular, we assume the following: 


Z(u)|T(-),0£» ~ Normal {y(u),a|(u)}; ueDsUDa, 


(14) 


where cKu) > 0, and 


y(u) =jU-hT(u;^J-f5(u;|); ueDsUDa, 


(15) 


is the unknown process of interest. In principal, one could easily adopt the generalized linear 
mixed effects model framework and replace the normal distribution in (fT4l) with the appropriate 
probability density function from the exponential class of distributions. For example, if Z(-) is 
count-valued than one might let Z(u) |y(u), be distributed as Poisson with the log link. 
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The unknown real value /i is interpreted as a eonstant “trend term.” Additionally, in (fTSl) we 
assume that a is an r-dimensional random veetor with mean zero and eovarianee matrix A = 
diag(Ai,..., A^). The speeifieation of the distribution of a, and assoeiated prior distributions 
for and A, are stated in Seetion 5. It is important to note that it is typieally straightforward 
to take an empirieal Bayesian approaeh by direetly estimating and A instead of plaeing prior 
distributions on these unknown quantities. 


The 5 nroeess rei 

oresents “fine-seale variabilitv.” We adont the models for 8 used in Wikle and 

Berliner ( 

2005 

) and 

Bradlev et al. 

(201.5b 

). That is, let ^ = {^j : j = l,...,ng)' eonsist of i.i.d. 


random variables with mean zero and varianee (7^, and let 


= (16) 

for any s G sueh that s is in the y-th areal unit inDg. Thus, = {\/\Bj\) fg.d(s;^)ds = ^j 

for Bj G Db- In general, (fT^ implies that the fine-seale variability term is eonstant within eaeh of 
the j = 1,..., ng areal units in Dg (with the respeetive value ^j). The speeifieation of the distribution 
of ^ and a prior for a'^ shall also be given in Seetion 5. 

3.2 CAGE for the Truncated Karhunen-Loeve Expansion 

It is not immediate that Proposition 2 (whieh motivated CAGE) holds for the proeess Y in (fTSl) . 
Thus, we provide an extension of Proposition 2 that develops the spatial aggregation error proper¬ 
ties of Y in (fTSl) . We formally state this result in Proposition 4. 

Proposition 4: Let f be any real-valued function with domain and be strictly greater than 
zero for each k= 1,..., r. Recall that a regionalization ofDg is given by Dc = {Cg : £ = 1, 
with C; nC^ = 0/or j ^ Ci ^ H C {1, ...,ng}, and Bh G Dg for £ = 1, ...,nc < ng. 

Define a generic point-level support {jc/ : j = l,...,nc}, such that Xj G Bj G Dg, where Bj C Cj 
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and j = Let = {Ysixj) : j = F^*'^ = (Fa(5^) : j= and Yc = 

{Ya{C) : A G Dc)'- Then the following statements hold for Y in 4771) .• 

fc) 

i. = ^{Cj; ^f)for j = 1, ...,no if and only if f{Y\ ’) = f{Yc) almost surely. 

iC') 

a. ^{Bj\ = 0(C;; ^fjfor j = 1, ...,nA, if and only iff{Y^g ’) = f{Yc) almost surely. 

iC') 

in. If ^^{xj) = ^{Cj\ ^ f for every Xj G Bj and j, then f{Yg j = fiYc) almost surely. 


Remark 7: For the process F in (fTTI) to have no spatial aggregation error on Dc we (again) require 
between scale homogeneity of the eigenfunctions. There are two key differences between Proposi¬ 
tions 2 and 4. The first difference is that Proposition 4 can be seen as an extension of Proposition 2 
from the multiscale K-L expansion in ® to the truncated process F in (fTTI) . The second difference 
is that Proposition 4 can be seen as a discretized version of Proposition 2. That is, Proposition 2 
allows Bj to be any subset of Aj, and Proposition 4 requires Bj to be defined on the (discrete) areal 
support Db. 


Remark 8: The choice to set r < oo is intimately related to the concept of spatial aggregation 
error. It is well known that predictors based on spatial basis functions with r-large disp lay more 


fine- level details than predictors based on spatial basis functions with r-small (ISteinL 


et ah, 


2013 


^Bradley 


2014af) . Thus, if r is chosen to be “too small” then predictions of F^ will have less variability 


over Ds (i.e., be more constant), and consequently the differences between F^ and Ya (or CAGE; 
see Proposition 3./) will be smaller than they should be. We strongly recommend performing an 
in-depth sensitivity analysis to choose r when using CAGE. To investigate the consequences of 
choosing r “too small” we provide a small sensitivity study in the Supplemental Materials. Addi¬ 
tionally, in the Supplemental Materials we provide a sensitivity analysis for the choice of r for the 
application in Section 5. 
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Similar to Proposition 2, we have that Proposition 4 guarantees that spatial aggregation error 
does not oceur for the spatial proeess in (fT5l) when a finite number of point-level eigenfunetions 
are eonstant over eaeh region in Da- This leads naturally to a definition of CAGE for the spatial 
proeess in (fTSl) : 


CAGE(A)=E 

DCAGE(C)=E 


f {0.(s) - 0.)yA{<|>,(s) - (^(A; , 1^' 

Ja |A| ' 

' ^ - 0(C; - 0(C; 


(17) 

(18) 


where A is a generie areal unit (i.e., A C Dfi, A = diag(A/ : i = l,...,r), C = UheH^h, H C 
{l,...,ne}, Bh E Db for eaeh h E H, and the expeetation is taken with respeet to the posterior 
distribution derived from (fT4l) and (fTSl) . Notiee that (fTTI) and (fTSl) are the truneated versions of 
CAGE and DC AGE in ® and ([7]), respeetively. In a similar manner a truneated version of Propo¬ 
sition 3 exists. We state and prove this result in Supplemental Materials. 


3.3 A Two-Stage Regionalization Algorithm 


The CAGE(A) measure allows us to evaluate whether or not the generie areal unit A has poor 
spatial aggregation properties. However, it is not immediately elear how it ean be used to speeify 
an optimal spatial support. We now deseribe the use of CAGE to explieitly obtain an optimal 
regionalization. Reeall that Db is the finest level aggregate support on whieh we wish to prediet. 
In general, our approaeh is to eonsider many different regionalizations (eombinations) of elements 
of Db and seleet from among them the support that produees the smallest average CAGE. By 
“regionalizations of Db” we mean a generie set Dc = {Ci : i — 1, ...,n^}, where CjHCi — % for 
j ^ i and for eaeh Cg — // C {1, ...,nB}, and Bg, E Db- 


A greedy seareh algorithm that seeks the minim um of the average CAGE (i.e ., CAGE(C£) /ng) 


poses a eonsiderable eomputational ehallenge (see 


Spielman and EoganI (12013h for related diseus- 
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sion). To address this computational issue we use a two stage search algorithm. In the first stage, 
a naive clustering algorithm is applied to each of the M samples of Yg from [Ye|Z], denoted 
Yg”^, for m = 1, For example, we could apply a /:-means algorithm to Y^”^ to define a set 
^(^)(yM) ^ : £ = 1,where is the £-th cluster returned by the /:-means algorithm. 

The superscript “(/:)” denotes the number of areal units in and we keep track of the depen¬ 
dence of the m-th replicate Yg^l In this article, we consider using the k-means algorithm. We set 
the input of the k-means algorithm to be the centroids of the areal units in Dg a nd Yg”l In the Su p- 


plemental Materials we also consider structural hierarchical clustering (SHC) (IMarslandl. 


20091) in 


place of k-means. The choice of clustering algorithm depends on the application. In settings where 
computation is of particular interest k-means is preferable over structural hierarchical clustering. 
However, structural hierarchical clustering allows one to incorporate neighborhood information to 
obtain contiguous areal units, which is a preferred regionalization in some applications. 

The first stage of our algorithm defines a collection of “candidate” spatial supports 


:k = gL, ...,gu-,m = 1, 


(19) 


Here, gi (gu) represents the smallest (largest) number of areal units one is willing to consider, and 
both gi and gu must be pre-specified. Notice that there are a total of M x [gu — gL + 1) spatial 
supports in which is considerably fewer than the total number of possible candidate spatial 
supports to chose from. 

In the second stage of the search algorithm we find the best (i.e., smallest average CAGE) 
subset of To do this, we compute 


= arg min 


£cage(c| 


i=\ 


( 20 ) 


where = {C^ '■ 7 = 1, and C for k 


It should be noted that D^, by 
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definition, is optimal since it is obtained by minimizing error. However, one might obtain a smaller 
value for the average CAGE by optimizing over a different set than Furthermore, one has to 
determine for their application whether or not it is appropriate to use CAGE or DCAGE in (l20l) : 
that is, in the case where the process is not interpretable on then one should replace CAGE in 
(l20l) with DCAGE. A step-by-step presentation of the regionalization procedure is provided in the 
Supplemental Materials. 

4 A Class of Multiscale Eigenfunctions 

Propositions 2 and 4 show that between scale differences in the eigenfunctions indicate that spatial 
aggregation error is present. Thus, the importance of the eigenfunctions for quantifying spatial 
aggregation error suggests that it should be parameterized. This will allow us to estimate eigen¬ 
functions, and hence, CAGE can be informed by the data. Below, we discuss the construction of 
what we call Obled-Creutin (0-C) eigenfunctions as a weighted combination of generic GBFs. We 
then discuss the properties of these basis functions. 


4.1 Obled-Creutin Eigenfunctions 

It has become common to express spatial random processes in terms of a b asis expansio n on ran- 


dom effects. As such, there are manv nossible choices for 

^asis functions ( 

Wikle, 

2010 

; Bradlev 

et ah. 

2014a 

). The insight provided by 

Obled and Creutin 

(1986 

) is that one can use any of these 


classes of point-level spatial basis functions to build an eigenfunction. We define an Obled-Creutin 
(0-C) eigenfunction as any real-valued function on Ds that takes the following form: 


(^)OC(s;F) = Y.Xj/i{s)Ftk; seD,,k=l,...,r, (21) 

i=l 
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where F is an r x r matrix with (i,fc)-th element given by the real value weight Fik^ and the r- 
dimensional veetor yf{-) = with —)■ M for i = l,...,r, eorresponds 

to the aforementioned GBF basis veetors. One ean organize the 0-C eigenfunetions into the r- 
dimensional veetor, F) = F),F))', whieh we eall an Obled-Creutin (0-C) 

veetor. 

It is not neeessarily true that Y (•; in (fTTI) leads to a multiscale truncated K-L expansion. 

In Proposition 5 below, we specify the condition such that 7(■; 0®*^) admits a multiscale truncated 
K-L expansion. 

Proposition 5: LetY j-; be the multiscale spatial process defined in 4731) . where Xj >0 

and > 0 for at least one j = l,...,r. Here, (•),..., V^r(-) are r real-valued functions with do¬ 
main Ds- Additionally, let F be an invertible r x r real—valued matrix. If F'WF = I then 
y j-; admits a multiscale truncated K-L expansion, where I is an r x r identity ma¬ 
trix and we define the (/, j)-th element of the rxr matrix W as Wij = ytj{s)ds. 


Remark 9: Proposition 5 is crucial for implementing the two-stage regionalization algorithm. 
That is, with a given GBF (i.e., radial basis functions, Fourier basis functions, wavelets, etc.) one 
can construct eigenfunctions, which can then be used within the two-stage regionalization algo- 
rith m from Section 3.3. There are many choices of GBFs available in the literature te.g.. Bradley 


et al. (2015at) l. and in Section 5 we use the local bisquare functions from ICressie and .lohannesson 


(120081). In the Supplemental Materials, we also consider using Wendland basis functions ( Wend- 


land. 


19981). 
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4.2 Specification of the O-C Weight Matrix, F 


We capitalize on the fact that the rxr matrix F is unknown. Estimating F will allow the data to 
inform the value of CAGE. However, Proposition 5 suggests that one needs to specify F with care; 
specifically, we require F'WF = I to ensure that y^(- ; 0®^) is a multiscale truncated K-E expan¬ 
sion. We achieve this by introducing a novel class of F matrices. This contribution is formally 
stated in Proposition 6. 


Proposition 6: For a given r-dimensional vector of basis functions iff let W be positive definite. Let 
G be an r X r real-valued orthogonal matrix. Then, 

F(G)=PwAw/"g, (22) 

satisfies F{G)'WF{G) = I, where is the Cholesky square root of the matrix 

Remark 10: For a given set of spatial basis functions { t///} we suggest verifying that W is positive 
definite. Then from (fT3l) . (|2T]) . and (l22l) one can write as 




^0C{-;F(G)}'a = V^(-)'F(G)a = V^(-)'PwAw/^Ga, 


(23) 


where a has mean-zero and rxr covariance matrix A. If a closed form expression for W is not 
available then numerical integration or direct Monte Carlo sampling can easily be applied to ap¬ 
proximate W. In the case of the latter, one can randomly generate n^, points {Sk'.l<.= \ ,...,n„} C 
using a uniform distribution on Ds, and approximate Wun with (l/nw)Lfc=i V^m(Sfc) 

In our Bayesian implementation given in Section 5, we use the following equivalent reparam¬ 


eterized expression of 




derived from the representation of in 




( 24 ) 
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where V^*(s)' = V^(s)'PwA^^'^^ for s G Ds, iff*(A)' = Xi Vi^Y^s for A G Da, and rj 

(= Ga) has mean zero and r x r covariance matrix Q = GAG^ Additionally, we assume that Q 
consists of random parameters that can be sampled. (For each application we undergo independent 
sensitivity analyses to select a prior distribution. For details behind the prior specification, and 
for related empirical results, see Supplemental Materials.) Then, it is straightforward to obtain 
samples of Q and tj, respectively, via a MCMC algorithm. Note that if a closed form expression 
for V^(s)^<is is not available then numerical integration or direct Monte Carlo sampling can 
easily be applied to obtain an approximation. In the case of the latter, one can randomly gener¬ 
ate points {sk : k = C A C Ds using a uniform distribution on A, and approximate 

■py Ia with (l/nw)Efeli In general, we have found that the value of needs to be 

large for these approximations to be reasonable (in Section 5 we set = 20,000). 

Additionally, one can obtain samples of the eigenfunction j-; F(Gf'"^)| to use within the 
expression of CAGE in ®. That is, denote the m-th replicate of Q with and let the cor¬ 
responding spectral decomposition be written as = Gf'”^A|j^Gt'"]^ Then, the corresponding 
m-th replicate of j-; F(G["*]) | is given by 


^°^|-;F(gH)| = i//*(.)'g['"]; m= 1,...,M. 


(25) 


We shall henceforth use the representation of Tj 

tion 0“|•;F(Gl’"l)| in 03. 




in (l24l) . and the 0-C eigenfunc- 


5 Application: Median Household Income from the American 
Community Survey 

We revisit the ACS 5-year period estimates of median household income for 2013 presented in Fig¬ 
ure 1. This data can be downloaded at http: //f actf inder2. census. gov/, This is an important 
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exam ple because there has been a growing interest in regionalizing data from ACS (' Spielman and 


Logan, 


2013 


20151). 


For this example, = 0, and = Da consists of the n = 3,109 counties in the continental US. 
Since US counties are the finest spatial resolution of the dataset in Figure 1, we set Dg = Da- Let 
[Z(-) |y (•)] be a normal probability density function with mean 7(■) and known variance Oz(-) >0, 
which are computed from margin of error estimates that are publicly available. Here, Z(-) is the 
log median household income, and we let 7(■) be distributed according to (fT5l) . 

Both a and ^ are assumed to be Gaussian, and we perform regionalization using V^), 


where ly(-) = : / = 


dCressie and Johannesson . 


. 1 5)' is a 75-dimensional vector consists of local bisquare functions 


2008b: 


{1-(I|S-C;||/W)2}2 if ||S-C^-|| <W 


wA^) = 


(26) 


0 


otherwise; s e Ds, 


with j = 1, ...,75 equally spaced knots Cj, and where w is 1.5 times the smallest distance between 
two different knots. The placement of knots was achieved using a space filling design t Nvchka 


1998b. We performed empirical studies that explore the relationship between r 


and Saltzman, 

and rffA (see discussion in Remark 8). These investigations suggest that r = 15 is appropriate for 
this example. (From our experience, our method is rather robust to the placement and number of 
knots, and the empirical results guiding this experience are provided in the Supplemental Materi¬ 
als.) We considered many different choices of prior distributions for the r x r cov ariance matrix 


0, and through independent sensitivity studies we found that the so-called MI prior (IBradlev et al.[ 


2015c 


appeared to be the most appropriate choice for this example (see the Supplemental Materi¬ 


als for more details). The k-means algorithm is used to define ^ in (fT9l) . and we let gi = 175 and 
gu = 195. Since the latent field is not interpretable on D^, we use DCAGE within the expression 
The variances of {e(A/) : / = 1, ...,n} are estimated a priori by ACS, and hence, are 


of D‘A’ in 
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(c) Optimal Support-Level DCAGE 



(d) State-Level DCAGE 




Figure 2: In (a), we present maps (for the eontiguous US) of predieted median household ineome 
(US dollars) defined on the optimal spatial support (i.e., F>^^) eonsisting of 185 areal units. Reeall, 
we eonsider areal units 175 through 195, and the value ehosen using DCAGE is 185. We superim¬ 
pose the state boundaries as a referenee to eompare to Figure 1(b). In (b) and (e), we present maps 
of the posterior standard deviations and DCAGE. In (d), we plot DCAGE by states. 

assumed known. 

In Figure 2(a) and 2(b), we present the predietions and eorresponding predietion error of me¬ 
dian household ineome on the optimal spatial support (and add state boundaries as a referenee). 
In Figure 2(b), the predietions appear fairly preeise with largest predietion error oeeurring in re¬ 
gions near Virginia, whieh have posterior standard deviation around 2,500 (whieh is roughly 5% 
of the mean median household ineome). The problems with spatial aggregation error indieated by 
Figures 1(a) and 1(b) deseribed in the Introduetion are no longer present in whieh eonsists of 
185 areal units. For example, eounties near Riehmond eonstitute a distinet region. Also, the state 
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of New York is divided into multiple distinct regions: areas near and in Manhattan, western New 
York, and upstate New York are all separated. However, it is worth noting that in Figure 2(c) the 
square root DCAGE values are comparatively larger around the state of Virginia. 

The DCAGE can also be used for uncertainty quantification. That is, state-level representatives 
may not be interested in the optimal regionalization produced by the two stage search algorithm, 
and instead, be interested in the median income over states. The DCAGE can be used to identify 
which states have poor spatial aggregation error properties. In Eigure 2(d), we plot DCAGE over 
states (i.e., treat states as fixed areal units), which has an average DCAGE of 0.24. This value is 
larger than the average DCAGE corresponding to the optimal solution, which is 0.19. Notice that 
the DCAGE corresponding to Virginia (and states near Virginia) are relatively high, while other 
states in the Midwest and West coast have comparatively smaller values of DCAGE. This would 
suggest that one should be concerned about assuming that statistics over Virginia can interpreted 
at lower spatial resolutions. 


6 Discussion 


The ecological fallacy and MAUP have become popular pedagogica ' 


ranhv and snatia 

statistics ( 

Robinson. 

1950; 

Onenshaw and Tavloi 

1979; 

Cressie. 

1993 

L Cressie 

and Wikle, 

2011; 

Baneriee et ah. 

201f 

5). However, very little has been done to characterize and 


tools for discussion in geog- 


mitigate these forms of spatial aggregation error from a statistical perspective. Thus, we provide a 
measure to formally characterize such error and a principled way to obtain an optimal (in terms of 
spatial aggregation error) regionalization defined over the generic continuous domain Re- 

gion alization has traditionally been solved using techniques outside the realm of statistics (Duou e 


et ah. 


2012 


Spielman and Eoganl 


2013 


Eolch and Spielmanl. 


2014 


Spielman and Eoganl. 


2015h . 


and our work offers a new perspective that respects the uncertainty of spatial random processes. 
Consequently, our methodology can significantly impact federal statistics, survey methodology. 
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geography, spatial statistics, and remote sensing/data acquisition settings. 

The heart of our methodology lies in the criterion for spatial aggregation error (CAGE), which 
we minimize to obtain our optimal regionalization. The methodological development of CAGE is 
intricate and involves a novel multiscale Karhunen-Eoeve (K-E) expansion. The introduction of a 
multiscale K-E expansion provides an approach to spatial COS that is not based on assumptions of 
between scale homogeneity. Eurthermore, the multiscale K-E expansion leads to a powerful tech¬ 
nical result that shows that any statistic does not suffer from spatial aggregation error as long as the 
multiscale eigenfunctions are homogeneous across scales. Thus, CAGE represents a measure of 
between scale homogeneity of eigenfunctions within a multiscale K-E expansion. There are many 
additional motivating features of CAGE, including connections to prediction error and across scale 
homogeneity of variances. 

To apply CAGE we need a parameterization of the multiscale eigenfunctions. This allows the 
eigenfunctions to be estimated, and hence, the CAGE can be informed by the data. Thus, we pro¬ 
vide a new class of Obled-Creutin tO-Cl eigenfunctions motivated bv the seminal naner of Obled 


and Creutin (Il986t) . The proposed class of 0-C eigenfunctions has broad applicability in the sense 


that any class of generating basis functions (GBE) can be used to build eigenfunctions. 

Einally, CAGE is used within an efficient two-stage regionalization algorithm. In the first stage 
of the algorithm (for a given number of areal units) a deterministic clustering algorithm is applied 
to each of the M samples from the posterior distribution of the latent process. This defines M 
spatial supports to select from. Then, in the second stage, the spatial support with the smallest 
(average) CAGE is chosen. This approach is extremely efficient, and accounts for the variability 
of the data by performing the search algorithm within the latent process space. 

An illustration of our algorithm was given using American Community Survey (ACS) 5-year 
period estimates of median household income. Comparisons of the optimal spatial support to the 
state-level ACS estimates indicate that the optimal regionalization preserves the county-level spa¬ 
tial information. Additionally, the size of this dataset is 3,109, and notably, the optimal spatial 
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support consists of just 185 areal units. The dramatic decrease of the dimensionality of the prob¬ 
lem has important implications for modeling very large spatial datasets. 

The applieation of CAGE to reduee the dimensionality of spatial data is just one of many ex- 
eiting avenues for future researeh. For example, the introduetion of spatially varying eovariates 
into the statistieal model will undoubtedly effeet the spatial aggregation error properties. Also, 
as previously mentioned, model seleetion eonsiderations, sueh as the number of basis funetions 
and elass of basis funetions, may effect the conclusions of the two-stage regionalization algorithm. 
The truncation of the multiscale K-L expansion is espeeially important from the point of view of 
regionalization, sinee fewer basis funetions lead to less variable predietions of the latent proeess, 
whieh then leads to fewer areal units produeed by the regionalization algorithm. Another inter¬ 
esting idea for future researeh would be to eonstruet a prior distribution for the regionalization by 
using the values of the CAGE to define prior weights. 

There are minor modifications to CAGE and the two-stage regionalization algorithm that would 
be reasonable to consider. For example. Proposition 2 shows that spatial aggregation error does not 
oeeur when point-level eigenfunetions are eonstant over eaeh region in the aggregate-level spatial 
support. Thus, we use the squared distanee between point-level and aggregate-level eigenfunetions 
to measure departures from the absenee of spatial aggregation error. However, other distanees be¬ 
sides the squared distanee might be used. This is similar to considering other forms of prediction 
error besides squared error. Also, there are a number of alternative seareh algorithms that one 
might consider. For example, one could use CA GE within a forward selection algorithm, or per¬ 


haps, one might use 


Spielman and LoganI (l2015[) ’s ACS regionalization (AReg) algorithm within 


the first stage of the two-stage algorithm. It would be diffieult to ineorporate AReg into the two- 
stage algorithm praotieally, sinee it is not computationally efficient for high-dimensional spatial 
datasets. The speeifications we use are eomputationally efficient and are shown to give favorable 
results. 
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I Introduction 


In this supplement to “Regionalization of Multiseale Spatial Proeesses using a Criterion for Spatial 
Aggregation Error,” by J.R. Bradley, C.K. Wikle, and S.H. Holan, we give additional insight to 
CAGE and the two-stage regionalization algorithm outside of what was presented in the main text. 
In partieular, we have applied the algorithm to another dataset, performed many different sensitivity 
analyses, and provided additional material that is meant to aid readers interested in implementing 
our proeedure. 

This supplement is organized as follows. In Seetion II, we provide guidanee on the implemen¬ 
tation of our algorithm ineluding: a summary of the statistieal model used in Seetion 5, details 
on prior distribution considerations, a step-by-step outline of estimation and the two-stage region¬ 
alization procedure, and additional discussion on model and regionalization specifications. Note, 
we use Roman numerals for section titles in this Supplement to distinguish from section titles in 
the main text. In Section III, we provide sensitivity analyses including: a comparison to a curren t 


state-of-art method for regionalization within the geography literature from 


Speilman et al. 


(12013h . 


a sensitivity analysis to the choice of Da, and a simulation study investigating the choice of the 
rank of the truncated multiscale K-E expansion. Next, in Section IV we provide a demonstration 
of the two-stage regionalization al gorithm to a dataset consisting of Mediterranean wind measure¬ 


ments (a subset of the data used in 


MilliffetaL 


(1201 IIB . which is used illustrate that the two-stage 


regionalization algorithm is flexible enough to handle multiscale spatial data. Einally, in Section V 
we provide the proofs to the technical results from the main-text. 


II Additional Details for Implementation 

Here, we give guidance on the implementation of our algorithm including: a summary of the 
statistical model used in Section 5 (Section Il.i), details on prior distribution considerations (Sec- 
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tion Il.ii), a step-by-step outline of the estimation and the two-stage regionalization procedure (Sec¬ 
tion Il.iii), and additional discussions on model and regionalization specification (Section Il.iv). 


Il.i Outline of the Statistical Model 


The statistical model introduced in Section 3.1 is summarized in Algorithm 1 below. We choose to 
describe this Bayesian hi erarchical model using the data, process, and parameter model terminol¬ 


ogy from 


Berlinen (Il996t) . 


Algorithm 1: Outline of the statistical model introduced in Section 3.1 


Data Model: Z(u)|/i,77,Q, ~ Normal{/r (u)'77 + 5(u;.§),a|(u)}; 

Process Model 1 : T] |Q ~ Gaussian (0, Q); 

Process Model 2: Gaussian ^0, ; 

Parameter Model 1 : ~ Normal ^0, ; 

Parameter Model 2 : cl ~ IG (a^,/3^); 

Parameter Model 3 : Q ~ [Q]; u G UD^. 


Here, the ng-dimensional random vector ^ = (i^i,...,> 0, > 0, /3^ > 0, and we let 

[Q] denote a probability density function for the unknown rxr covariance matrix Q. We consider 
many different choices for [Q], and provide these details in Section Il.ii. The value of is chosen 
to be large so that the prior distribution on /i is interpreted to be vague, and similarly, we set 
oc^ = ^ so that the prior distribution on is flat. 
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Il.ii Prior Distributions to Consider 


As ISorbye and Rud (120141) discuss, the prior distribution (and the assoeiated hyperparameters) on 
the r X r eovarianee matrix Q affeets posterior inferenee. As sueh, we eonsider several different 
ehoiees for priors on eovarianee matriees. In partieular, we eonsider three different prior distribu¬ 
tions. The first prior distribution we consider is the eonjugate inverse Wishart distribution. This 
is a fairly common choice because it allows for direct sampling of the full-eonditional dis tribution 


corre sponding to Q, however, in high-dimensions this prior is known to perform poorly (IHodges 


2013b . 


The seeond prior distribution we eonsider is from 


Bradley et al. 

(201 4 J 

) and 

Bradley et al. 


(l2015cb . where it is assumed that 


Q = ^[R«‘.e^+{Qi(l-A)Q*}Rs'] 


( 2 . 0 ) 


where is the best positive approximate (IHighaml. 


1988b of a square real-valued matrix M, 


<7^ > 0 is unknown, the ngx r matrix 'Pg = (V^(5)': B G Db)', = Qb^b is the QR deeompo- 

sition, and A is the ng x ng adjaeeney matrix eorresponding to Z)g. Notiee that (jb]) ineorporates 
spatial information, but is not spatially refereneed. That is, this prior for Q is motivated by spee- 
ifying eov('Pg?7) so that it is “close” to the eovarianee from an ICAR model on Dg, where 'Pg 
is spatially refereneed but T] is not. An inverse gamma prior i s placed on where t he hyperpa- 


rameters are ehosen b ased on the suggest i ons in S eetion 3.2 of 


Bradley et al. 


(12014bb and 


Bradley et al. 


Sorbye and Ruel (l2014b . Following 


(l2015cb . we refer to this prior speeifieation as the “MI” 


prior distribution due to a eonneetion to the Moran’s I statistie. 


The third prior d istribution we consider is the Givens angle prior (lYang and Berger , 


Bradley et al. 


1994 


2015bb . where the speetral deeomposition is written as Q = EqA^Pq, and the rxr 
diagonal matrix Aq has diagonal entries set equal to the eigenvalues of (©. The parameter is 


assumed to follow a flat inverse gamma distribution (i.e., with shape and seale set equal to 1). The 
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r X r orthogonal matrix Pq is decomposed into a Givens rotator product, 


Pq = (Ol ,2 X Oi ,3 X • • • X Oi^r) X (02,3 X • • • X O 2 ,,-) X • • • X Or-i,r, 

where Ot j is a r x r identity matrix with the (i,i)-th and {j,j)-th element replaced by cos{Oij) 
and the (/, 7 )-th (( 7 ,/)-th) element replaced by —sin(0;j) (sin( 0 ,-^ 7 )). Here, 6ij G [—:;r/2, ;r/2] is 
unknown, and let the shifted and rescaled 6ij be denoted as Qj = 1/2 + Oijjn. Then, it is assumed 
that 

=a + bxgij{PQ); i<j = l,...,r, (2.0) 

where logit(^,- 7 ) = log{^/,y/(l - a,/? G M, and gijCPq) represents the Givens angle 

of Pq. Finally, a vague Gaussian prior is placed on {a,by (i.e., Gaussian with mean zero and 
variance 1000). For all of our analyses we considered all three prior distributions. These sensitivity 
analyses suggested that the MI prior lead to the best predictive performance for the application in 
Section 5, and the inverse Wishart prior led to the best predictive performance in Section V. 

Il.iii Outline: Estimation and Implementation of Regionalization 

In this section we give a brief outline of the two-stage regionalization algorithm. It should be 
acknowledged that, for any given application, minor modifications to these steps may be needed. 

1. Define the spatial support Dg, which represents the finest resolution one is willing to predict 
on. If = 0 we suggest setting Dg = Da, which is the finest resolution information that 
is available. When 7 ^ 0 then one has the freedom to choose any spatial support for Dg, 
however, one should be mindful of the size and spatial coverage of the locations within Df. 
Thus, for illustration, when 7 ^ 0 we suggest setting Dg to a fine resolution grid. 

2. Obtain M MCMC replicates of Yg = {Ya{B) : B G Dg)', using the statistical model in Algo¬ 
rithm 1. Specifically, let rj represent the m-th replicate of T] and | represent the m-th 
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replicate of (§. Then, the m-th replicate of Yg can be computed as 


where the ngx r matrix = (l//‘*(u)': u G Dg)'■ The Bayesian procedure can easily im¬ 
plemented using a Metropolis with in Gibbs sampling algorithm. 

3. Use a naive clustering algorithm to obtain ^ in (19). We consider two clustering algorithms 
to define namely, the fc-means algorithm, and structural hierarchical clustering. In gen¬ 
eral, the fc-means algorithm takes on as it’s argument an ng x / real-valued matrix J, and 
returns a clustering of the rows of J. Let L be a ng x J matrix with the y-th row equaling 
the centroid of the y-th areal unit in Dg. Then, we let / = J -f 1 and set J = [L, Yg^^]. The 
structural hierarchical clustering approach takes on two arguments J = [L, Yg”^] and the 
adjacency matrix corresponding to Dg. 

4. Choose the spatial support from that minimizes CAGE. That is, compute according 
to (20). If Y can not be interpreted on substitute CAGE with DCAGE. 

5. Produce maps of the values in the sets {Ya{C^p) : G {var(}( 4 (C''^^|Z)): g 

and {CAGE(C^P) : C^p g (or {DCAGE(C''^^) : C^p g when appropriate). This 
allows one to visualize the process and its corresponding prediction and spatial aggregation 
errors. 

Il.iv Model and Regionalization Algorithm Specifications 

To implement the two-stage regionalization algorithm, we need to specify: the number and place¬ 
ment of knots that define the r-dimensional GBE and the lower and upper bounds on the number 
of areal units used within the two-stage regionalization algorithm (i.e., gg and gu). We now pro¬ 
vide discussion on to make these choices in practice. 
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Specification of Knots: The choice of knots and r is important for preserving the appropriate fine- 
scale features of 7^. If the fine-scale features of Ys are ignored then the two-stage regionalization 
algorithm may produce too coarse of a regionalization (see simulation study in Section IV.iii). 
However, the number of areal units produced by the two-stage regionalization algorithm appears 
to be robust to r “too large.” Recall the number of areal units in is denoted with This 
interaction between the number of optimal areal units and r suggests an approach for selecting the 
rank r, which we outline into the following steps: 


(1) Consider a fixed range of values for r (i.e., r = rj/). 

(2) For each r = r^, use the algorithm outlined in Il.iv to find an optimal regionalization 
and n'^. There will be a different value of for each each r = r^, 


(3) Plot r versus n^. 

(4) Choose the value of r to be the point in which does not change dramatically as r in¬ 
creases. 


We follow the suggestion of iRuppert et al. 


(2003 


, chap. 13, pp. 255-260) and apply a space filling 


design algorithm to a set of randomly selected points {cj : j = 1,..., r*}, wher e we set r* = 600 > r . 


The space-filling design can be determined using the FUNFITS function in R (INvchka et al. 


1998 1. 


Then, we choose r according to steps 1—3 above. For the applications in Section 5 and Section V 
we found that, respectively, r = 15 and r = 200 are appropriate. 


Specification of gi and gu’. The widest range of values that we can consider for regionalization is 
gL = '^ and gu = ng —I. To specify less extreme choices for gL = 2 and gu = ng — I we consider 
running a simplified version of the two-stage regionalization algorithm, and use the results of the 
“simplified two-stage regionalization algorithm” to inform a tighter range between gu and gu- 
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In particular, we first run the two-stage regionalization algorithm (outlined in Seetion Il.iii) with 
M = I, gi = 2, gu = n — and use the means algorithm. Then, we ehoose gi and gu to be a 
tight range eentered around found using this simplified two-stage regionalization algorithm. 


Ill Simulations, Sensitivity Analyses, Comparisons, and Techni¬ 
cal Clarifications 


Here, we provide many different side-studies ineluding: a simulation study to eompare the two- 
stage regionalization algori thm to a eurrent state-of-the-art alternative in the geography literature. 


Spielman and LoganI (1201 51) ’s ACS regionalization (AReg) algorithm (Seetion Ill.i); a small sensi¬ 


tivity analysis on the ehoiee of Da (Seetion Ill.ii); and a simulation study investigating the ehoiee 
of the rank of the spatial basis funetion expansion (Seetion Ill.iii). 


Ill.i Simulation Study: A Comparison to 


Speilman et al. 


( 2013 ) 


In this seetion, we establish that our approach performs regionalization extremely well relative to 
the AReg algorithm available in the geography literature. To do this, we generate synthetie data 
based on a subset of the ACS 5-year period (from 2009 to 2013) estimates of the pereentage of 
households below the poverty threshold. We generate the spatial field. 


Z{A)=Ya{A)+£{A)-, AeDa, 


(3.0) 


where Da is the set of 351 census traets surrounding the eity of Austin (TX). Let {Z(A)} represent 
the perturbed version of the logit transformed pereent below the poverty level ACS survey estimate 
(denoted by { 1 ( 4 (A)}). (Notiee that we use the symmeterizing logit transformation, where, for a 
given pereentage p, logit( 7 ?) = pl(\-p).) The set {e(A) : A G Da} eonsists of independent normal 
random variables with mean-zero and known varianee. The published varianees for percent below 
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the poverty level are transformed to the logit seale using the delta method (lOehlert . 


1992h . and 


used as the known variances of {e(A)}. In practice, the ACS estimates (i.e., {1^} for this example) 
are publicly available and are, hence, observed. Nevertheless, for the purposes of this simulation 
study we will act as if the ACS estimates are an unobserved spatial field to be estimated from Z. 
To obtain we model this data using the mixed effects model in Algorithm 1, where V^(-) = 
: j — 1. ....42)' is a 42-dimensional vector consists of local bisouare functions t Cressie and 


Johannesson, 


20081): 


I {1-(||S-C;||/W)2}2 if||s-Cj||<W 
Xj/jis) = < (3.0) 

0 otherwise; s G 


with j = 1,..., 42 equally spaced knots Cj, and w is 1.5 times the smallest distance between two dif¬ 
ferent knots. Note, that we are not restricted to using local bisquare functions, since our modeling 
framework is general enough to allow for any desired GBF. For computational convenience, we use 
the k-means algorithm to define ^ in (19), and let gL = '^ and gu = 100. The latent process in ® is 
not defined on Ds, and thus, we shall use DCAGE within the expression of in (20). Addition¬ 
ally, we denote the output of AReg with : k = 1,..., and compute it using 

software made available at https: //github. com/geoss/ACS_Regionalization/blob/master/README.md. 

The goal of this simulation study is to compare the error properties of and This is 

done using the following metrics: 


ReMSPE(ZA) 


ReCAGE(ZA) 


£?= r { Ya (Af") - (Af")} ' 

op /■ N 2 


AReg . ^ 

tA, cAf") 

\^\ 


Id 
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where /(•) is the indieator funetion. Here, ReMSPE stands for “relative mean squared prediction 
error” and ReCAGE stands for “relative spatial aggregation error,” respectively. Values of ReMSPE 
that are larger (smaller) than 1.0 indicate that prediction on has smaller (larger) MSPE than 
when predicting on Thus, values of ReMSPE that are larger (smaller) than 1.0 indicate that 

the two-stage algorithm (AReg) leads to better (worse) predictive performance. Eikewise, values 
of ReCAGE that are larger than 1.0 indicate that the two-stage algorithm is preferable in terms of 
spatial aggregation error. 


(a) Histogram of ReMSPE 

25 n. 



0.9 1 1.1 1.2 1.3 1.4 1.5 


ReMSPE 


(b) Histogram of ReCAGE 

25-| 



0 5 10 15 20 25 

ReCAGE 


Eigure 3: In (a) and (b), we present histograms of ReMSPE and ReCAGE from taken over the 100 
replicates of Z defined in (21). The red line indicates the value of 1 in both panels. A value of 
ReMSPE and ReCAGE greater than 1.0 indicates that the two-stage regionalization algorithm is 
preferable over AReg. 


We simulate 100 replicates of Z in ®, and compute ReMSPE and ReCAGE for each of the 
100 replicates. Eor both metrics our proposed algorithm consistently outperforms AReg. In fact, in 
each of the 100 replications of Z we obtain an ReMSPE > 1.0, and a ReCAGE > 1.0, where ReM¬ 
SPE ranges from 1.0112 to 1.3979 and ReCAGE ranges from 5.8408 to 23.1620, respectively (see 
Eigure 1 for a histogram over the 100 replications of Z). It is somewhat expected that ReCAGE 
suggests that the two-stage regionalization algorithm is preferable over AReg because from Propo¬ 
sition 3, CAGE is directly related to the squared difference between the lower spatial resolution 
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process and the aggregate-level estimator. However, it is rather interesting that ReMSPE suggests 
that the two stage algorithm is also preferrable in terms of squared predietion error, sinee AReg is 
motivated by redueing sampling error. This may be due to the faet that AReg does not take into 
aeeount survey error (i.e., {e(A)}), while the two-stage regionalization algorithm aeeounts for this 
error by performing its seareh in latent spaee. 

Ill.ii Sensitivity to Da 

Notiee that the two-stage seareh algorithm takes on Ds and Da (the spatial domains of interest) as 
an input. Thus, one might be interested in the sensitivity of our approaeh to the spatial domain of 
interest. For example, in Figure 2(a) we plot the optimal areal units (i.e., found in Seetion 5, 
over California, Oregon, Nevada, and Arizona. Now, suppose we let Da eonsist of the 126 eounties 
in California, Oregon, Nevada, and Arizona, and we re-run the two stage seareh algorithm on this 
restrieted domain (i.e.. Da no longer eonsists of all eounties in the mainland of US, but eonsists 
only of eounties in California, Oregon, Nevada, and Arizona). The D°(? found under this restrietion 
is given in Figure 2(b). 

There are 12 areal units in without restrieting Da, and 11 when one restriets Da- Upon 
eomparison of Figures 2(a) to 2(b) we see that the general pattern of the two-stage seareh algorithm 
is robust to this ehange in Da, however, the final answer does ehange. We note that sinee the 
initialization of the k-means algorithm is random, the eandidate set of areal units are not neeessarily 
the same eaeh time one runs the two-stage seareh algorithm. 

Ill.iii Simulation Study: Selection of the Rank of the Truncated Multiscale 
K-L Expansion 

In this seetion, we use simulation to investigate the impaets of misspeeifying the rank of the trun- 
eated multiseale K-F expansion. In partieular, we ehoose a simulation model with r = 100 random 
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(a) equal to all US counties (b) equal to CA, OR, NV, AR only 



Figure 4: In (a), we plot the optimal areal units (i.e., D^^), found in Section 5, over the state of 
California. In (b), we plot the found by restricting Da to consist only of counties in California, 
Oregon, Nevada, and Arizona. Each distinct color identifies a different areal unit, and the relative 
difference between each color is arbitrary. The state boundaries are superimposed as a reference. 

effects, and we perform regionalization with r misspecified and r correctly specified. The region¬ 
alization with r correctly specified is treated as the “correct” regionalization, which we compare 
to. 

Let the latent process of interest be generated as follows: 

T,(s) = n + + 5(s; 4); s G (3.-1) 

where Ds = {s = (5i,5'2)^: ■s'iA 2 = [0.05,0.1,0.15,..., 1] x [0.05,0.1,0.15,..., 1]},recall Ys{s;(j>^^) 
is defined in (11), and let be based generated from 100 equally spaced (over Ds) local bisquare 
basis functions. The corresponding dataset is generated as follows: 

Z,(s) =Ysis) + £sisy,seDf CDs 

Za{s) =YA{A)+eA{Ay,AeDA, (3.-1) 

where we randomly select 50% of the observations from Ds to define D^, and Da consists of the 
10 X 10 grid cells that cover [0,1] x [0,1]. We let ej(-) be a mean zero white-noise process with 
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Point-Level Data 


Process 



Figure 5: Example simulated data and process. These maps are produced using (I3.-1I) and (13.01) . 
The top left panel contains simulated data on (with 50% of the field being covered). The top 
right panel contains the simulated process on D^. The bottom left panel contains the aggregate data 
process (i.e., Z^), which has complete spatial coverage over Da- The bottom right panel displays 
Ja. 

variance = 0.1820 (so that the signal-to-noise ratio (=5) is large). Likewise, {eA(^) : ^ G Da} 
consists of i.i.d. independent mean zero random variables with variance 0.1820, and is independent 
of the spatial random process ej(-). An example of the data and the process is given in Figure 3. 

Consider performing regionalization using the outline in Section Il.iii, to the data in Figure 3 
with r = 9,100, and 256. For illustration let Dg = Da, and set gi = 2 and gu = 99 (the largest 
possible range). Here, r = 9 represents the case where r is too small, r = 100 represents the case 
where r is correct, and r = 256 represents the case when r is too large. When r is too small we 
obtain fewer areal units (6 areal units) than when r is correct (13 areal units); however, the optimal 
regionalization algorithm is robust to the case where r is too large, which produced 15 areal units. 
This conforms to intuition as it is well known that predictors based on spatial basis functions with 
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r-large display more fine-level details than pred i ctors b ased on spatial basis functions with r-small 


(IBradley et ah. 


2011 


Stein. 


2013 


Bradley et al 


2014ah . Thus, one would expect that if r is chosen 


to be “too small” then predictions of will have less variability over (i.e., be more constant), 
and consequently lead to coarser regionalizations. 

These conclusions are similar over multiple replications; in Figure 4 we provide histograms 
of rff? obtained from the two-stage regionalization algorithm over 50 independent replications of 
{Zs} and {Za}. Notice, however, that the variability associated with r too large is much higher 
than when r is too small and when r is correct. The p-value of a sign test comparing when 
r = 9 (r = 256), to n'^ when r = 100 is 0.0494 (0.5716), which suggests that when r is too small 
(too large) we obtain coarser (similar) results than when r is correct. 

The fact that there is no significant change in the number of areal units when r is too large also 
conforms to intuition; since there are enough spatial random effects to capture fine-scale behavior, 
and the remaining random effects are negligible. This interaction between the number of optimal 
areal units and r suggest an approach for choosing r (i.e.. Steps 1—3 in Section II.iv). For the 
ACS application in Section 5, we consider r = 25,50,75,100,125, and 150. Likewise for the 
Mediterranean wind example we consider r = 25,50,75,100,125, and 150. In Figure 5, we plot 

versus r (i.e.. Step 3 from Section Il.iv). Here, we see that for the applications in Section 5 and 
Section V we found that, respectively, r = 15 and r = 200 are appropriate. 

Ill.iv Technical Clarifications: Positive Definiteness of the Multiscale K-L Ex¬ 


pansion 


A covariance function cov{ 1 ( 5 ( 8 ), Tj(u)} is positive definite if dCressie 


1991 p. 68), 


L L ^i^jCov{y5(s/),Ti(S;-)} > 0 

i=i;=i 


(3.-1) 
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Histogram of optimal number areal units (r =; 9) 

0.81 -^^^^^— 



optimal number areal units 


Histogram of optimal number areal units (r = 100) 

0.51 -^^^^^- 



optimal number areal units 


Histogram of optimal number areal units (r = 256) 



optimal number areal units 


Figure 6: Histograms of over 50 independent replieations of {Z^} and {Z^}. The value of r 
used to fit Algorithm 1 is indicated in the title of the panel. 

for any finite number of spatial locations {s, : / = 1, and any set of real numbers {bi : i = 
1,m}. That is, the covariance function, associated with the spatial random process Y^, is positive 
definite if a weighted average of covariances implied by any set {Tj(s,) : i = l,...,m} has non¬ 
negative variance, where {bi : i = 1, are the generic weights. The validity of the covariance 
of Ti in (11) follows immediately from the definition of positive definiteness, and the quadratic 
form of 

cov 
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n°P versus r (Section 5) 


n°P versus r (Section iV) 
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Figure 7: The plot of versus r as described in Section Il.iv. In the left panel we plot versus 
r for the ACS example presented in Section 5, and in the right panel we plot versus r for the 
wind example in Section IV. The values of r considered in the ACS example in Section 5 were 25, 
40, 50, 75, 100, and 150. The values of r considered in the wind example in Section IV were 50, 
75, 100, 150, 200, and 250. 


where A is defined below Equation (15) of the main text, 

and 

That is, let b = {bi,...,bm)', and notice that 

171 171 

£ £ bibjcov {y,(s,), } = cov = b''p('”)A'p('")'b > 0, 

and hence, ® holds for the covariance associated with in (11). In a similar manner, one can 
prove the validity of the covariance function of T in (1) using Proposition l.ii. 
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IV Application: Mediterranean Surface Winds 


A critical component of the interfaee between the atmosphere and the upper oeean oeeurs due to 
the transfer of momentum and the exehange of heat and fresh water, whieh is manifested through 
surfaee winds from the atmosphere. Due to a laek of direet measurements of surfaee wind over the 
oeean, wind field estimates over sueh regions were historieally based on a blend between meeha- 
nistic models of the atmosphere and a relatively sparse global network of wind observations from 
buoys and ships of opportunity. The praetieal spatial resolution of these so-called “analysis” winds 
is limited to fairly large spatial and temporal seales of variability, yet they are reported on fairly 
high-resolution grids. The advent of spaee-bome seatterometer instruments in the 1990s provided 
the first high-volume, high-resolution in spaee, wind estimates over the oceans. Although these 
seatterometer winds have higher spatial resolution (effeetively “point” seale), they a re inc omplete 


in space and time, necessitating an optimal blending approach (e.g.. 


liff et al. (1201 lb . and 


Wikle et al. 


Wikle et al. 


mom. Mil- 


(120131) give reviews of recent statistical approaches to generate 


spatially and temporally complete ocean wind fields. 

As mentioned above, the weather center analysis winds do not contain spatial infonnation 


commensurate with the spatial support in which they are estimated (e.g., see 


Milliff et al 


(1201 ll) 


for discussion). That is, the kinetic energy spectrum of the winds does not contain realistic variation 
at small spatial scales. The support given by the additional (and incomplete) seatterometer wind 
estimates is relatively much smaller. To date, there have been no attempts to consider an optimal 
spatial support for statistical wind predictions given these types of data. 

In the example presented here, we consider ocean surface wind data from two sources over the 
Mediterranean Sea. In particular, we consider the north-south wind component for analysis winds 
from the European Center for Medium range Weather Forecasting (ECMWF) and satellite wind 
observations from the QuikSCAT seatterometer; this is a subset of the data used in the study by 


Milliff et al. 


(1201 ih . We assume that the high resolution (25-km) seatterometer wind observations 
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Figure 8: Wind observations from 2 February 2005 at 12:00 UTC (Universal Coordinated Time), 
(a) North-south (v) component of the wind from the ECMWF-analysis winds on a 0.5° x 0.5° 
grid, (b) North-south wind component from the high resolution (25km), but spatially intermittent, 
QuickSCAT scatterometer wind retrievals. 


are effectively “point” support (relative to the analysis winds). Thus, these data are recorded on 
both Ds C and Da- Here, Ds ranges from 30° to 48° north latitude, and -19° to 42° east 
longitude, and Da consists of a 0.5° x 0.5° resolution grid on Ds. In total. Da consists of 4,551 
areal units and Ds consists of 6,916 observations for the time of interest, resulting in a dataset of 
11,467 spatial observations. Figure [8] shows these data for a 6-hour window centered on 12:00 
UTC (Universal Coordinate Time) for 2 February, 2005. 

In this application, we let Dg be a half-degree grid. We consider the model in Algorithm 1, 
where l/A is a multiresolution bisquare basis vector consisting of local bisquare functions in ®. 
We chose r = 200 knots using a space-filling design and the plot in Figure 5 (see Section Il.iii). 
We consider both structural hierarchical clustering and k-means to define ^ in (19) with gi = 280 
and gu = 380; note that we these choices of gi = 280 and gu = 380 were guided by the approach 
discussed in Section Il.iii using the k-means algorithm with initial choices of gi = 2 and gu = 600. 
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We also considered an equivalent analysis using the Wendland GBFs with fc-means clustering. 


Here, the Wendland basis functions (IWendland . 


19981) are defined as 


v^r(s)= 


(1-JXs))^(35JXs)^ + 1H(s)+ 3)/3 ifO< < 1 
0 otherwise; s E Ds 


(4.0) 


where j = 1,...,200, (ij(s) = ||s —c*||/w, we choose w = 1.5 times the smallest distance between 
two different knots, and {Cy} consists of the same 200 knot specifications used in the bisquare 
basis functions. Additionally, since th e latent field i s inter pretable on Ds, we use CAGE within the 


expression of in (20). Following 


MilliffetaL 


(1201 lb . the variances of e(u) are set equal to 1 


when u G Ds, and set equal to 10 when u E Da- 

The results of the CAGE analysis of the posterior wind predictions is given in Figure HI The 
top row of this figure shows that when using the standard 0.5° resolution support, there is a notice¬ 
able high CAGE “crescent” in the south central portion of the region. This would suggest that one 
should be concerned about assuming that statistics on the wind field over this region can be inter¬ 
preted at the point level. Note that the optimal support regions with k-means and bisquare GBFs 
(the second row of H]) are much larger than the Dg level shown in the first row, but the predic¬ 
tions look qualitatively similar to the half-degree predictions, although with more smoothing and 
the corresponding reduction in root prediction error associated with the relatively large optimal 
aggregation regions. The optimal aggregation seems to pick up realistic meteorological features. 
For example, notice the homogeneous region centered on Corsica and Sardina, which corresponds 
to a region of more intense southerly winds off of the mainland (so- c alled “Mistral winds”) that 


are important in forcing the ocean circulation (e.g., see 


MilliffetaL 


(1201 lb ). Perhaps more im¬ 


portantly, although the higher CAGE crescent is still present, it is noticeably reduced in intensity 
relative to the Dg support. The Wendlend GBF predictions (third row) are similar to the bisquare 
predictions, but with generally larger regions and with higher CAGE values that are shifted north- 
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ward. Finally, the last row of Figure [9] shows the bisquare results with the struetural hierarehieal 
elustering method. These are similar to the bisquare fc-means results, but one notices more spatial 
detail in the predictions. 


There is a striking amount of dimension reduction that results from the CAGE analysis. That 


is, values of are considerably smaller than the number of observations, 11,467. We have that 


Hfi = 323 when using the bisquare GBFs and k-means, n^i =315 when using the Wendland GBFs 
and k-means, and = 327 when using the bisquare GBFs and SHC. This suggests that optimal 
aggregation, such as the results presented in Figure 7, may be a viable alternative approach for 
dimension reduction. 

We note that there is quite a large amount of shrinkage in these wind predictions relative to 
the data, which is not surprising given the uncertainty in the winds and th e fact that no temporal 


information is being considered here. As discussed in 


Wikle et al. 


(12013h . one can gain signifi¬ 


cant prediction efficiencies if temporal dynamic information is included in the model for winds. 
Such an analysis is beyond the scope of this simple illustration, but the CAGE-based selection of 
prediction support could, in principle, be utilized in that framework. 


V Technical Proofs 

In Section V.i, we provide the proofs to Propositions 1—6. In addition to these proofs, we also 
provide results alluded to, but not explicitly stated in the main text (Section V.ii). 

V.i Proof of Propositions 1—6 


Proof of Proposition 1: The assumptions of Pro position 1 allo w us to apply the K-L de 


composition of {T(s) : s G from 


KarhunenI (Il947h . That is, from 


KarhunenI (Il947h we have 
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Half Degree Grid Predictions 


Half Degree Grid Root Squared Prediction Error 


CAGE on Hait Degree Grid 




-6 

-8 
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-10 



Optimal Support-Level Root Prediction Error 



CAGE for Optimal Support 





0.4 


-6 

-8 


0.3 


-10 
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Optimal Support-Level Predictions (Wendland) 




CAGE for Optimal Support (Wendland) 



Optimal Support-Level Predictions (SHC) 

r 


Optimal Support-Level Root Prediction Error (SHC) 




CAGE for Optimal Support (SHC) 



Figure 9: CAGE-based posterior summaries of the predieted north-south wind eomponents based 
on the analysis and seatterometer observations from 2 February 2005 at 12:00 UTC. The first 
eolumn displays the posterior mean; the seeond eolumn displays the posterior standard deviations; 
and the third eolumn eontains the ealeulated CAGE. In the first row the values (i.e., posterior 
mean, posterior root predietion error, and CAGE) are all defined on a half degree grid. In the 
seeond row values are defined on the optimal spatial support found using /:-means and the bisquare 
GBFs. In the third row values are defined on the optimal spatial support found using fc-means 
and the Wendland GBFs. In the fourth row values are defined on the optimal spatial support 
using struetural hierarehieal elustering (SHC) and bisquare GBFs. Note that the eolorbar for the 
predietions differ from the eolorbar used in Figure 6. 


20 














































that for se Ds 

oo 

YA{Bh) = Y,H^)0Cj. (5.0) 

j=i 

where the eigenfunetions {<^>;(s) : j = 1,2,...} have domain Ds and satisfies, 

/ <t>ji.^)<Pk{^)ds = 5jk, (5.0) 

J Ds 

where djk is the Kroneeker delta funetion. Additionally, the random variables in the set {aj : j = 
1,2,...} are uneorrelated with varianees {Xj : j = 1,2,...}, and the coefficients {ccj : j = 1,1,...} 
can be found by projecting ys(-) onto the eigenfunctions. That is. 


aj= / y,(s)(^);(s)Js, 


ID, 


(5.0) 


for eac h j. Also, these eigenfunctions are solutions to the Fredholm integral equation ('e.g. JPapoulis 

lod)), 

C{s,\i)<pj{s)ds = Xj<pj{\i)-, \xeDs,j =1,2,..., (5.0) 


IDs 


where, from the statement of Proposition 1, C(s,u) is a valid covariance function for each s,u G Z)^. 
The statement that 

CO 

Ya{A) = Y,M^)(^d ( 5 - 0 ) 


i=l 


in L?-{Q.) for A C Ds, is equivalent to saying that 


UA)^e{ iYA{A)-'£(^A,j{A) 


a. 


(5.0) 


i=\ 


converges to zero as n goes to infinity. Note that in ®, the expectation is taken with respect to 
(f2, 1^). Expanding ® we have. 


C„(A)=£{1'.,(A)2}+£<^ E'i’AyW 


V /=1 


a,- 1 y - 2 e\Ya{A)[Y, ^aj{A) aj 

0=1 


(5.0) 
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The first term of the right-hand side of ® can be written as 


E{YA{Af}=E!^^-^I^Jj,{s)Y,(u)dsdu\ 

The second term of the right-hand side of ® can be written as 


^ 1 ( H ^A,j(A)CCj 
,!=1 


) V/=i 


^11 0a,/(A) 0aj(A) 

U'=if=i 

W^ee/J, ^s,ii^)^s,j{yi)0Ciajdsdu 
1 " " /■ /■ 

uaHH ^s,i{s)^sj{n)E{aiaj)dsdu 
1^1 ,-=i^=i^aJa 


since recall from the K-L decomposition that a, and «/ are uncorrelated with variances A,- and A/, 
respectively. Finally, the third term of the right-hand side of db]) can be written as 


E 
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Since a, is found by projecting onto the eigenfunctions. From ® we have that 




0=1 


= E 


^F.(w)0/(w) JwJsJuj 
j^Ys{u)Ys{w)(j)i{w)dwdsdu 


s) / £’{yj(u)yj(w)}0,(w)JwJs<iu 


wLLp-'^ 

^ L ^ '^) (w)^^wJs Ju. 


From the Fredholm integral equation in ® we have 


i = J^C{u,yf)^i{\v)dyfdsdu 


^ ^ E (^),,,-(s)(^,,/(u) A/JsJu. 


Substituting ®, ®, and db]) into © gives 


j^C'(s,u) 


( 5 .- 7 ) 


Upon taking th e limit as n goes to infinity on both sides of ®, it follows from Mercer’s theorem 


jMercei . IQOol) that 


lim^„(A) =0, 

n—T’oo 


( 5 .- 7 ) 


for each A C note that Mercer’s theorem shows uniform convergence at the point-level, allow¬ 
ing one to pass the limit through the integral. This proves the result. 
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The proof of l.ii follows a similar logic to ®. That is, note that 


£ ^aM) MB)^i - cov {Ta (A), Ta (5)} 


/=1 




Upon taking th e limit as n goes to infinity on both sides of ®, it follows from Mercer’s theorem 


jMercer . 19091) 


19091) that Proposition l.ii holds. 


Proof of Proposition 2: First, we prove the following statement: If for 

('A') 

j — 1, ...,nA and for any positive integer k, then - Ya almost surely. Then the continuous 

mapping theorem is applied to get /(Y^’^^) = /(Ya) almost surely. 

fA) 

We proceed using a proof by contradiction. Assume that Y^ is not almost surely equal to Ya. 
Then, for at least one X; and A/, there exists a 7 > 0 such that 


P (| T ,( x ,)- TA ( A ,)|> r )> 0 . 


(5.-8) 


However, we have from Chebychev’s inequality 


p(|y,(x,-)-TA(A,-)|>r)< 


{T,(x,-)-TA(A/)r 


yZ 


(5.-8) 


Assume that 0yt(x;) = ^A,kiAj) for j = l,...,nA and every positive integer k. Then, upon adding 
and subtracting ^ki.'^i) within ® we have: 


P(|T,(x,) -Ta(A, 0| >y)<^E\ y,(x,-) - £ U^i)a,+ £ ^A,k{Ai)ak-YA{A, 


T 


k=\ 


k=\ 


n 


—E I y,(x/) - Y, hi'f^doCk > + I Y <^A,k{Ai)ak - YA{Ai 


r 

2 

+ ^E 

yl 


k=\ 

n 


yZ 


k=\ 


Ysi^i} Y^ ^fe(Xi)0ffc / s Y^ ^A,k{Ai)CC]^ YA(Ai 


k=\ 


k=\ 
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It follows from iKarhunenI (Il947h that the first term on the right-hand-side of ® eonverges to zero. 
Likewise, from Proposition 1 the seeond term on the right-hand-side of dH) eonverges to zero as n 
goes to infinity. Note that sinee P(|ys(x,) — YA{Ai) \ > 7 ) does not depend on n we have that, 


p(|y,(x,-)-y^(A0|>7) 


< lim -^E 

n—5.00 yZ 


- 52 h{^i)cCk \ < Y. ^A,k{Ai)ak-YA{Ai) 

k=i ) U=i 


Thus, we are left to find the expression of the limit in ® . Note, 


n 


^2 ^ki^i)C)Ck r \ ^2 ^A,k(Ai)OCi( Ya{Ai) 
k=l j [a:=1 

Yys(xi)ak^k(s)ds 


±El£r^(^,)r,(u)du 


~ I E E (l>ki'^i)(l>ji^)cCjCCkds 


IC=lj = 
n 


+ { E ^k{'f^i)^kYs{s)ds 


For the term in (jh]) notiee from (jh]) and © we have that 




£y,(x/)^ y,(u)(^)fe(u)Ju ^k{s)ds 

]T\ f i f E{Ysixi)Ysiu)}(l)j{u)du (l)kis)ds 

^ / E / C{xi,u)(j)k{u)du (j)k{s)ds 
|A| JAi ^“1 Jd, 


If" 

TJT Y^k{s)h{Xi)?ikdS. 

|A| 
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The terms in ® and ® ean be written as 


- r'" {ly.My.(u)du^ = -Tfi C(x,u)<;u}, 

\ \ " f 1 If" 

~ u\^ ^ E E / <l>ki's^i)<l>jis)0Cjakds V = - — / £ <^fe(s)0^(x/)AfeJs. 
1^1 I A:=l I 1^1 


For the term in ® notiee from ® and ® we have that 


7 ^ e \ ^ [ 4(x,')a^y,(s)Js i = ^ ^ 0fe(s)4(x,')A^Js. 


Thus, it follows that 


Y,{xi) - (pki-fi-dcck \ < (l)A,k{Ai)ak-YA{Ai) 


k=l 


k=l 


2 f " 

= m / E <^>fe(X;)0/t(s)Afe-C(x/,s)Js, 


whieh, again by Mereer’s theorem, eonverges to 0 as n goes to infinity. Thus, from ® we have 
that 

F(|T,(x,)-TA(A,)|>r)= 0 , 

whieh contradiets ® . One can prove forward implication of Proposition 2.ii in a similar manner. 

To prove the reverse statement of Proposition 2.i, suppose that = /(Ya) almost surely 

for any measurable real-valued function /. Thus, the functions fi{h) = bi for i = l,...,nA and 
b = {bi : / = 1, ...,nAy G W^, imply that 


T,(xO = Ta(A/), 


(5.-12) 


almost surely. Multiplying both sides by cty we have 

T^(x;)ttj = }a(A,-)(Z^- 
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almost surely. Substituting ® into the equation above gives, 


Ys{xi) [ Ys{s)(j)j{s)ds = -^ [ [ Ys{u)Ys{s)(j)j{s)dsdu. 
Jd, \Ai\jAiJD, 


Taking the expectation on both sides we have 


[ C{xi,s)^j{s)ds =[ [ C{u,s)^j{s)dsdu, 

JDc \-^i\ JAiJDs 


and then from ® we have 


Dividing by A/ 


(^j{xi)Xj = ^j{u)duXj. 


This proves the result. One can prove the reverse statement of Proposition 2M in a similar manner. 
By the condition in Proposition l.iii, we have that for a given 


UBj) = 


^ 0 fc(s)(is = 


\Bj\ Jbj 


Bj\ Jb 


^k{Aj)ds 


= (l>k(A 


1 




j\ ’’Bj 


Ids = ^k{Aj 


It follows from Proposition 2.ii that Proposition 2.iii holds. 


Proof of Proposition 3: We now prove the equalities listed in Equations (8), (9), and (10) 
of Proposition 3. We start with Equation (8). Notice that for a given s G Ds, A G Da, 
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{h}. 


{y,(s)-y^(A)}^|{ 0 ,},{A,} 

= E <!y.(s)-£0fc(s)a4 KM,{4} 


k=\ 


+ E 


k=l 




E2E 

+ 2E 

+ 2E 


^ki^)(^k ^A,k{‘^)C)Ck r I{ 4 }, {4} 

k=\ ) 

j^^A,k{A)ak-YA{A) \ km ,{ 4 } 


k=\ 


y^s) - £ h{s)ak j> <{ £ (l>kis)cck- ^ (j)A,kiA)ak [ K<^>/t},{4} 

fe=l J [fe=l A:=l J 


y^s) - £ (^)yt(s)a^ \ <{ £ (j)A,kiA)ak-YAiA) i KM,{ 4 } 

k=\ j 4=1 J 


£ (l)kis)ak- Y. ^A,k{A)ak \ < Y <l>A,ki^)0Ck-YA{A) \ K<^>fe},{Afe} 

.k=i k=i J 1^=1 


Through an application of Mercer’s theorem we have that the sum of the cross -produet terms in 


Karhunen ( 


1947 h 


dH), dH), and dll) eonverge to zero as n goes to infinity. Similarly, it follows from 
that dll) goes to zero as n goes to infinity, and from Proposition 1 that dH) goes to zero as n goes to 
infinity. Thus, 


f=l 


Then, upon taking the expeetation with respeet to {(pk}-, {4}|Z we have the desired result. 
To prove Equation (9) recall from Mercer’s theorem and Proposition l.zz that. 


var{yKs)} = Y^kisf^j 


k=l 


var{yA(A)} = Y (pAA^f^j- 

k=l 
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Expanding (9) and substituting ® we have 


CAGE(A) = E 


= E 


= E 


= E 


= E 


W ' 

L"., - 21"., ((.;(s)04,y(A)Aj 


| 2 l| 


ds+ ^ 0A,fc(A)^Aj|Z 


k=i 


\A\ 

L 7 =M^nj 

\A\ 

‘ var{y,(s)} 

I |A| 


Js - 2 £ + £ (j)A,k{Af^Z 

k=l k=l 

oc 

Js- £0A,,t(^)%iz 

k=\ 


(is —var{}A(A)} |Z 


; A CDs. 


This proves (9). 

We now prove Equation (10). Erom (8) we have for any A C Ds, 


CAGE(A) = 


{y,(s)-yAO)}- 

|A| 


■ds\Z 


(5.-19) 
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Adding and subtracting Ya, 


CAGE(A)= E 


= E 


y,(s)-yA(A)+yA(A)-yA(A) 


y,(s)-yA(A) 


|A| 

2 


-JslZ 


+ E 


-ds\7j 


Ya{A)-Ya{A) 


|A| 


-<is|Z 


+ 2E 


= E 


y,(s)-yA(A) yA(A)-yA(A) 


y,(s)-yA(A) 


|A| 

2 


■JslZ 


-JslZ 


+ E 


yA(A)-yA(A) Js|z 


IE 


= E 


yA(A)-yA(A) |z 


-1 

1 

NJ 

_1 



/ -1—1- —dsYL 

Ja |A| ' 

-E 

|yA(A)-yA(A)| |z 


This proves Equation (11). 

Proof of Proposition 4: The fine-seale variation term 5 in (16) ean be written as 

5 (u;(|) =h(u)'(|; ugD^UDa, 


where 



(/(u G 5) : 5 G Ds)' ifuGD, 
(M; 5 eZ)gy ifuGDA, 
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and /(■) is the indicator function. Then, from Equation (15) we have that for a given (j>^ and a. 


Yf^ = Atl„^ + 4>Fa + HF| 

Yc = Illnc + ^cCC+iic^, 

where the nc x r matrices = (0^(x^)': j = l,...,nc)' and 4>c = : j = 

(C) 

and the nc x ng matrices H) = (h(xj)': j = 1, ...,nc)^ and He = (h(Cj)': j = 1, Notice 

(C) 

that for the values of {x^} and {C^} given in the statement of Proposition 5, we have H) = He = 
(the nc X nc identity matrix), and thus, 

Yf^=Atl„^+4»Pa + ^ 

Yc = Atlnc+^ca + l- 

(C'^ I—I 

The condition for the forward implication of Proposition A.i is that 4*^ = 4>c; thus, from ® we 

have that 

yP = At 1.C + a +1 = Yc. (5.-27) 

When applying any real-valued measurable / to both sides of db]), we obtain that f{YP) = /(Ye) 

almost surely. One can prove forward implication of Proposition A.ii in a similar manner. 

(C) 

To prove the reverse statement of Proposition A.i, suppose that f{Y\ j = fiYc) almost surely 
for any real-valued function /. Thus, the functions //(b) = bi for i = 1, ...,nA and b = {bj : j = 
1 , ...,nA)' e imply that 

Yf^=Yc, (5.-27) 

almost surely. From (jb]) and © we see that 

4>Fa = 4>ca, (5.-27) 

almost surely. Multiply both sides of ® by a', and take the expectation with respect to T|0^, A to 
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obtain 


4>Fa = 4>cA. (5.-27) 

Provided that Xj > 0 for all j, we can take the inverse of A on both sides of ® so that = <E>c, 
which is the desired result. One can prove the reverse statement of Proposition 4.ii in a similar 
manner. 

By the condition in Proposition we have that for a given and a, 

^,(x,)'a = ^(C,; j=\,...,nc. (5.-27) 

Integrating ® with respect to xj we have 

<l>,ya = <l>{Cf, 0J'a; j=l,...,nc. 


Since Xj > 0 for all j, this leads to the condition for the forward implication of Proposition 4.ii, 
and thus, it follows that Proposition 4.iii holds. 

Proof of Proposition 5: From Equation (1) we see that for y(-;0®^) to be a multiscale 
truncated K-L expansion, we only need to sho w that is a tru ncated K-L expansion. 


Obled and CreutinI ( 19861) . 


Many of the following equations can be found in[; 

To show that Ys{-;(j>y^) is a truncated K-L expansion, we need to establish three items: the 
eigenvalues must be nonnegative with at least one eigenvalue strictly positive; the Fredholm inte¬ 
gral equations must hold; and the eigenvectors must be orthonormal. Notice that 


cov 




OC/ 


Wii^)FikOCk r ^ ^ Wqi^FqpCCp 

^k=li=l J {q=lp=l 

E 4 < £ ¥i{s)Fik \ < £ ¥qi^)Fqk 
k=l [/=! J U=1 
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Substituting ® into the Fredholm integral equation we have, for = 1,r, 



EEE 

i=lk=lq=l ) l^m=l 


ds= (Or, 


,?=1 


qp 


(5.-29) 


where {cOyt} represents the eigenvalues of Ys{-;^^^). Distributing the sums and integral through 
® , we obtain 


r 


q=l 




¥q{^)Fc 


qp 


Matehing terms in ®, we have 


(5.-29) 


FqkhFikWimFmp = (OpFqp-,q = 1, ...,r. (5.-29) 

i=l k=l m=l 

In matrix form, ® beeomes. 


FAF'WF = Fa, (5.-29) 

where A = diag(Ayt) and Q. = diag((Oyt)- The assumption that F'WF = I and ® implies that the 
Fredholm-integral equation holds provided that 

FA = Fa. (5.-29) 

Sinee, F is invertible we have that ® verifies that the eigenvalues of are nonnegative 

with A = a (and at least one eigenvalue is strietly positive), and that the Fredholm integral equa¬ 
tions for Ys{-',^f^) hold. The orthogonality of (j>^^ holds by assumption sinee 

/ (^)pc(s; F)<^)pC(s; F)Js= ^ L / v^fc(s)v/p(s)rfs 

k=\p=l 
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which results in the relation, 


F'WF = I. 


This eompletes the proof. 

Proof of Proposition 6: Let W = PwAwPw be the speetral deeomposition of W. It follows 


that the Cholesky square root of W and W ^ is given by PwA^^ and PwA^^'^^, respeetively. It 


follows immediately that G'(PwA^^'^^)'WPwA.^^^^G = I. 

V.ii Additional Results 

In the main-text, three results were discussed, but not formally stated. Thus, in this seetion we state 
and prove these results. In partieular, at the end of Remark 6, we mentioned that the CAGE identi¬ 
ties in Proposition 3 also hold for DCAGE; this extension of Proposition 3 is referred to as Result 
1. Also, at the end of Seetion 3.2 we mention that a version of Proposition 3 exists for CAGE in 
(17) and DCAGE in (18); these two extensions are referred to as Result 2 and Result 3, respeetively. 

Result 1: Assume that the conditions of Proposition 1 hold. Assume that the stochastic process 
Z : X —)■ M A generated based on any generic probability space (Q,^, 3^') such that the 
conditional probability density function ofY («) |Z exists for each u G DsUDa, where Z is defined 
in Remark 2. Then, DCAGE in (7) has the following alternative expressions: 


DCAGE iC) = E 


y {YA{Bh)-YA{C)f 
2-1 \r\ 


DCAGE iC) = E 



y var{TA(^/;)} 

h \c\ 


var{yA(C)}|Z 


DCAGE iC) = E 
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where C = // C {1, ...,ng}, and G Ogfor each hEH. 


Proof of Result 1: We now prove the equalities listed in Equations dH), dS), and dS) of Propo¬ 
sition 3. We start with Equation dH). Notiee that for a given G Db, C = UheHBh, H C{l,...,«s}, 
{(^^(•)},and {Xk], 


{YA{Bh)-YA{C)Y\{h}Ah] 


= E 


+ E 


YAiBh) - ^ ^A,k{^h)CCk > 


k=\ 


^A,k{Bl'i)CHk ^A,kiC^OCj^ / 


k=l 


k=l 


. 1 Z7 


+ 2E 

+ 2E 

+ 2E 


^MC)(^k-YA{C)\ \m,{h} 


k=l 


\ 


YA{Bh) ^A,k{Bh)^k r > ^A,k{Bh)(^k ^A,k{C)(Xj^ / |-[A^} 

k=l J U=1 k=l J 


k=l 


Ya(B,) - £ MBh)0Ck'^ ^t.MC)0Ck-YA{C) \\{h}dh} 


^A,k{^h)CCk- 52 ^A,k{^)CCk [■ 52 ^A,k{^)CCk-YA{C) > {{(pkjA^k} 

k=l k=l j [fe=l J 


Through an application of Mercer’s theorem we have that the sum of the cross-product terms in 
dH), dH), and dH) converge to zero as n goes to infinity. Similarly, it follows from Proposition 1 that 
dH) and dH) go to zero as n goes to infinity. Thus, 


{YA{Bh)-YA{C)Y\{h]Ah} = £ {^AABh)-^Aj{C)flj, (5.-33) 

f=i 
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Then, upon taking the expeetation with respeet to {0^}, {Afc}|Z we have the desired result. 
To prove Equation ® reeall from Proposition \.ii that, 

oo 

\dx{YA{Bh)] = ^ (^A,k{Bhf^j 
k=i 

oo 

var{yA(C)} = 

k=l 

Expanding ® and substituting ® we have 


CAGE(C) = E 




heH 


\C\ 


f 'Li=\^A,j{Bh)^Xj 2Y,j^i^A,j{Bh)^A,j{C)Xj ^ \2'i \rj 

-+ L 

heH k=\ 

Y^L7=1<^)A,;(5/0%- 

^ \ 1^1 ^\^^A,k{C) A/+ Y ^A,kiC) 2lj\X 

JieH 1^1 k=\ k=\ 

£ 11 

heH 1^1 k=\ 


= E 


E^^^%!^-var{r,(C)}|Z 


heH 


C 


A CD, 


This proves ® . 

We now prove Equation ® . Erom ® we have. 


CAGE(C)= E 


y {YA(BH)-YA{C)}^ ^^ 


heH 


C 


(5.-39) 
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Adding and subtracting Ya, 


CAGE(C)= E 


= E 


UB,,)- Ya (C) + Ya(C) - Ya (C) 
--- -1^ 

heH 


c\ 

2 


^ Ya(Ba)-Ya(C) 

Z- - 


heH 


C 


+ E 


Ya{C)-Ya{C) 

-^IZ 


heH 


C 


+ 2E 


= E 


YA{Bh)-YA{C)\\YA{C)-YA{C) 

-^iz 

heH 


ICI 


YA{Bh)-YA{C) 

Z- --^IZ 


heH 


C 


+ E 


Ya{C)-Ya{C)\ |Z 


2E 


yA(c)-yA(c) |z 



{Ya{B,,)-Ya(C)Y 


'l 2 

E 

heH 1 

-E 

yA(c)-yA(c) |z 





This proves Equation ® . 

Result 2: EorZ defined in (14) and Y{-; (j>fi defined in (13), we have that CAGE in (17) has the 
following alternative expressions: 


CAGE (A) 
CAGE (A) 


= E 


= E 


{y,(s; ^s)-Ya{a-^_^ 


Ja |A| 

f var{y,(5; 

Ia |A| 


ds-war{Y a{A; ^J}|Z 


CAGE{A)= E 


1 

(N 

1 

4^ 



- 1—1 - —ds\Z 

-E 

[Ya{A)-Ya{A- |Z 


where A is a generic areal unit (i.e., A C Dfi, and Ya{A) = E {yA(A) |Z}. 
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Proof of Result 2: We now prove the equalities listed in Equations ®® , and ® of Propo¬ 
sition 5. We start with Equation ® . Notiee that for a given s G Ds, A G Da, OL, 0^, and A, 

<I>s)-Ya{A; - (/►(A; <j>^)}. 

Taking the expeetation with respeet to OC|0^, A we have 

|Eii:[{n(s; 

(5.-45) 

Then, upon taking the expectation of ® with respect to 0^, A|Z and integrating s over A, we obtain 
Equation ®. 

To prove Equation ® notice that 


var{y,(s; = ^,(s)'A^,(s) 
var{yA(A; 0J} = ^(A; 0J'A0(A; 


Expanding ® and substituting ® we have 


CAGE(A) = E 
= E 

= E 

= E 

= E 




|A| 


A |A| 

j 0.(s)'A0,(s) -20,(s)^A^(A; 

[ 0.(s)'A0,(s) 

Ja |A| 

I 0.(s)'A0,(s) 


ds\Z 


-Js + ^(A; 0J'A0(A; 0J|Z 


|A| 

var{y,(s; 

|A| 


Js-20(A; 0J'A0(A; 0J + ^(A; 0J'A0(A; 0J|Z 
Js-0(A; ^,)'A0(A; 0J|Z 


Js-var{yA(A; 0J}|Z 


■, A CDs. 


38 











This proves ® . 

We now prove Equation ® . From ® we have for any A C Ds, 


CAGE(A)= E 


{Us-. 

W 


(5.-51) 


Adding and subtracting Ya, 


CAGE(A) = E 


= E 


y,(s; ^^)-Ya{A)+Ya{A)-Ya{M <I>,) 


|A| 


-ds\Z 


n(s; 0,)-D4(A) 

w 


-dslZ 


+ E 


Ya{A)-Ya{A-, 


|A| 


+ 2E 


= E 


-2E 


= E 


y,(s; 0 J-yA(A)}{y^(A)-yA(A; 


IAI 


■ds\Z 


n(s; ^,)-fA(A) 


-JslZ 


-JslZ 


+ E 


Ya{A)-Ya{A-, Js|Z 


Ya{A)-Ya{A-, |Z 


y,(s; <I>,)~Ya{A) 


-JslZ 


Ya{A)-Ya{A-, |Z 


This proves Equation ® . 


Result 3: Eor Z defined in (14) and Y (•; defined in (13), we have that DCAGE in (18) has the 


39 





























following alternative expressions: 


DC AGE {C) = eIy, \Z 

[heH 1^1 

DCAGE{C)= ^J) _var(y^(C; 0,))|Z) 

heH I 

DCAGE{C)= e[ £ \z\ -E \{Ya{C) - Ya{C-, (I>,))^\z\, 

[heH \E\ J '■ '' 

where C = Ufj^fiBh, H C {1, ...,ng}, and G Dgfor each h E H. 

Proof of Result 3: In the proof of Result 2, replaee the integral with sums, and replaee 0^(s) and 
7^(8;with and YA{Bh,(^,), respeetively. 
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